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abstract 


A method for including the solution of the transfer equation in 
a standard Henyey type hydrodynamic code has been developed. This 
modified Henyey method has been used in an implicit hydrodynamic 
code to compute deep envelope models of a classical Cepheid with a 
period of 12^ including radiative transfer effects in the optically 
thin zones. 

There are two secondary features on the light curve of the 
model, a shoulder during rising light and a distinct bump during 
falling light. It is shown that the shoulder during rising light is 
caused by a deep envelope pressure wave and that the bump during 
falling light may be due to an atmospheric oscillation. It is shown 
that the atmospheric oscillation mechanism is consistent with the 
Hertzsprung sequence and the period-luminosity relation. 

The structure of each hydrodynamic model was used as a snapshot 
of the temperature and pressure structure of the atmospheric layers. 
After using line blocking factors to account for the effect of 
spectral lines on the spectral energy distributions computed for the 
models, broad band UBVRI colors were calculated. The light and color 
curves of the models reproduce the observed amplitude and asymmetry 
of Cepheids in this period range. In addition, the color-T e £f 



relations derived from the models were found to agree with those 
derived independently* It was found that the colors of the 
equilibrium model are best reproduced if the intensity mean of the 
magnitudes, <B > - < V is used to compute mean colors „ It 

was also found that loops in the (U-B)-(B-V) diagram are probably 
due to the dependence of the continuous opacity on the electron 
pressure. 

Line profiles were then computed using the moving atmospheres 
from the hydrodynamic models . It was found that the velocity 
gradients in the atmosphere are not responsible for the large 
microturbulent velocities observed in Cepheids but may be responsible 
for the occurrence of supersonic microturbulence. The total observed 
microturbulence was found to be consistent with the linear sum of 
the classical microturbulence and that caused by the velocity gradients. 
It was also found that the splitting of the cores of the strong lines 
is due to shock induced temperature inversions in the line forming 
region. 

The adopted light, color, and velocity curves were used to 
study three methods frequently used to determine the mean radii 
of Cepheids. It was found that an accuracy of 10% is possible only 
if high quality observations are used. 
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NOTATION 


The choice of notation always presents a problem. One would like 
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These two criteria are often incompatible. In the following standard 
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meaning, the context in which it is used is sufficient to remove any 
ambiguity i The following is a partial list of the symbols used. 
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E radiation energy density 
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CHAPTER I 


INTRODUCTION 


A. Historical review 

In 1784, John Goodricke discovered the variability of § Cephei. 
However not until 126 years after he published his findings was the 
importance of Cepheids as distance indicators realized. Although it 
was known that Cepheids differed from other variables in that their 
velocity curves are mirror images of their light curves (Belopsky 
1895), it was not until 1912 that Leavitt discovered the period- 
luminosity relation (Pickering 1912). During Leavitt’s study of 
variable stars in the Small Magellanic Cloud, she noticed that the 
brighter stars have longer periods. Since all stars in the Small 
Magellanic Cloud are nearly the same distance from the earth, she 
concluded that the correlation applied equally well to the absolute 
magnitudes . 

Unfortunately, the early attempts by Hertzsprung (1913) and 
Shapley (1918) to find the zero point of the period- luminosity 
relation neglected the effects of interstellar absorption and were 
based on poor data. The resulting error of iTs almost exactly 
compensated for the error introduced by using globular cluster 

Cepheids to find the absolute magnitudes of RR Lyrae stars. When 

* 

the absolute magnitudes of RR Lyrae stars were found independently, 
the results supported Shapley' s zero point. It was not until 1952 
when Baade was unable to find RR Lyrae stars in M31- -stars that he 
should have been able to see — that the error was discovered (Baade 


1 



2 

1956). The importance of Cepheids as distance indicators was 
demonstrated by : this correction which increased both the distance 
and time scales of the universe by a factor of two. 

The discovery of the period-luminosity relation increased the 
theoreticians' interest in the nature of the Cepheid mechanism. 

While it was originally thought that Cepheids were binaries, 
difficulties reconciling this hypothesis with observations suggested 
that the light and velocity variations were due to radial oscillations . 
Although Ritter had suggested radial pulsations in 1879, and Shapley 
had introduced the idea into the astronomical literature in 1914, it 
was not until 1918 that a thorough study of linear, adiabatic pulsa- 
tions was published by Eddington . In this work, he derived the 
period -mean density relation, P V p 88 Q, a constant, and showed that 
the sign of the second order terms dropped from the linearized equations 
indicated that these terms were responsible for the observed asymmetry 
of the light curves. He was unable to determine the nature of the 
driving mechanism but suggested that an increase of the energy genera- 
tion at minimum radius would produce the desired effect, A more 
troubling problem was the cause of the 90® phase lag. Linear, adiabatic 
theory predicts that maximum light should be in phase with minimum 
radius, not with maximum expansion velocity as observed. If the light 
and velocity vary sinusoidally, this phase shift is 90° „ 

Partly because of these difficulties, geometric theories to explain 
the light and velocity variations were not abandoned. However, in 1926, 
Baade proposed a test of the pulsation theory. If the star is pulsating, 
the observations at two phases can be used in 


the relations 





3 

(1-1) 

(1-2) 


Since he had no knowledge of the color-temperature relation 
and the necessary velocity measurements had not been made, he was 
unable to apply this method . The first successful radius determination 
using this method was made by Becker in 1940, but the results were not 
very accurate, Wesselink (1946) improved Baade's method by selecting 
phases of equal color. If equal color was assumed to imply equal 
temperature, no temperature calibration was needed in equation (1-1). 
The radius determinations made using Baade's, and later, Wesselink's 
methods firmly established the validity of the pulsation hypothesis. 

The nature of the driving mechanism continued to be a problem. 

Four possibilities had been suggested: 

1. € mechanism: During contraction, the energy generation 

increases. Since heat is added to the star when it is 
hottest, thermal energy is converted into mechanical 
energy (Eddington, 1918) . 

2. . K mechanism: During contraction, the opacity in some region 

increases. Since less heat is lost at maximum temperature, 
the stability of the star is decreased (Zhevakin, 1953). 

3. y mechanism: An ionization zone will remain cool during 

compression since some energy will be used to ionize the gas. 
The gas will absorb heat when hottest, destabilizing the 
star (Cox, Cox, Olsen, King, Eilers, 1966). 



4. r mechanism: At minimum radius, the increased curvature 

of the outer stellar layers traps radiation. Heat is 
added at maximum temperature, and the stability of the 
star is reduced (Baker, 1967). 

The radius calculations of Becker (1940) indicated that the 
relative radius variations in Cepheids were too small for the r 
mechanism to be significant. In 1950, Epstein showed that the radius 
variations in the stellar core were so small that the e mechanism 
must be negligible. Linear, nonadiabatic calculations of Baker and 
Kippenhahn (1962, 1965) and by Cox (1963) demonstrated the effective- 
ness of both the K and y mechanisms in the Hell ionization region. 

B. Recent work 

There were other problems that could not be investigated with 
the linear approximation. The 'nonlinear calculations of Christy, J. 
Cox, A. Cox, and King, among others (see King and Cox 1968, and 
Christy 1969, 1970 for references) have been used to study these 
problems. These models show that the present theories are adequate 
to describe the gross features of the pulsation such as the approxi- 
mate light and velocity amplitudes, the phase lag between radius and 
temperature changes, and the shapes of the light and velocity curves. 
Other questions remain unanswered. Stobie's (1969c) calculations 
show that there is a line in the period- luminosity plane separating 
stars pulsating in the first harmonic from those pulsating in the 
fundamental. In general, stars with a period less the 7^ should be 
pulsating in the first harmonic; those with a period over 7^, in the 
fundamental. On the other hand, Fernie (1968) has shown that the 
scatter in the empirical period -radius relation is reduced if some 



Cepheids are treated as overtone pulsators. The stars that Fernie 
suggests are overtone pulsators are not confined to the period range 
predicted by Stobie. The discrepancy may be due to errors in the 
radius determinations inherent in the Wesselink method. These 
errors will be investigated in Chapter VI. 

The nature and cause' of the phase lag are also problems. 

Castor (1968) has suggested that the phase lag is caused by the 
hydrogen ionization region moving through mass as the star pulsates. 
Using results of linear theory, he concludes that the heat capacity 
of this region delays light maximum. He states that the rate at 
which the ionization region sweeps through mass should be in phase 
with the luminosity. Christy (1968) has studied this problem with 
his nonlinear calculations and finds that the phase shift through the 
hydrogen ionization region is only 30°. The remainder of the shift 
he attributes to skewing of the light curve by nonlinear effects. 
King, Cox, Eilers, and Davey (1973), on the other hand, find a 90° 
phase shift in their linear calculations and conclude that the phase 
lag is a nonadiabatic effect and not a nonlinear effect. The phase 
lag will be investigated further in Chapter III. 

Probably the most worrisome problem is the Cepheid mass 
discrepancy. Stellar evolution calculations show that the mass 
of a star near the Cepheid instability strip has a mass given by 

A- 

M ev ~ ex P t 2 « 3 (X + 3Z) ] 

(Iben and Tuggle, 1972a). Masses of pulsating stars, on the other 
hand, can be found in one of two ways. The first uses the period- 
mean density relation, P Vp = Q. Linear calculations can be used 



to find Q (Epstein 1950; Cogan 1970; Cox, King, and Stellingwerf 1972) 
and, if the radius can be found, the relationship yields the pulsation 


mass, Mq. It should be noted that is sensitive to errors in the 

3 

radius, since Mq ~ R . The second method uses the phase of the second 
bump on the light curves of Cepheids with periods between 7 d and 15^. 
Christy (1968) has proposed that this bump is the result of a pressure 
wave which travels into the star, is reflected from the core, and 
appears at the surface as a second bump on the next pulsation cycle. 
The time it takes the pressure wave to travel to the core and back 
is a measure of the stellar radius, while the period of the star de- 
pends on both the mass and radius. It is possible, therefore, to 
find the mass from the phase of the bump, M^. Unfortunately, and 
are typically half of 

Iben and Tuggle (1972a) have shown that the descrepancy between 

M and can be removed by increasing the zero point of the period- 
ev Q 

luminosity relation by 0^2 or adjusting the re l a tion 

slightly. Fricke, Stobie, and Strittmatter (1971, 1972), however, 
have shown that the discrepancy cannot be removed if is also 
considered. Since van Genderen (1970) has suggested that more than 
one mechanism produces bumps on Cepheid light curves, M| may be un- 
reliable. This point is examined further in Chapter III. 

More complete reviews of the literature have been given by 
Rosseland (1949), Ledoux and Walraven (1958), King and Cox (1968), 
and Fernie (1969). 


C. Hydrodynamic Cepheid atmospheres 

One of the difficulties with the models discussed above is that 
they do not adequately represent the Cepheid atmosphere. Since the 
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atmosphere is the only part of the star actually observed, the coarse 
zoning and use of the diffusion approximation in the optically thin 
zones of the models makes the comparison of the theoretical and observa- 
tional results difficult. Recently Keller and Mutschlecner (1970, 1971) 
and Bendt and Davis (1971) have attempted to remove some of these 
difficulties by including the effects of radiative transfer. Their 
emphasis, however, is still on the envelope. In Chapters IV and V, an 
attempt will be made to answer some of the following questions raised 
by the observations: 

1. How should light and color variations be averaged to best- 
represent the equilibrium state of the star? For example, 
there is a systematic difference between 


2 . 


3 . 


1 

V 



(B-V)dt 


and 


-2.5 log 


X io ~ B/2 ' 5 dt 

/V v/2 - 5 dt 


Why are weak lines asymmetric near phases of maximum velocity 
while strong lines are often asymmetric near phases of minimum 
velocity? For example. Bell and Rodgers (1964) find that the 
X4508A line of Fell is asymmetric in j3 Dor near the phase of 
maximum radius. 

Why are the cores of the strongest lines sometimes split? 
Grenfell and Wallerstein (1969) and Wallerstein (1972) observe 
splitting in the core of that indicate velocity differences 
of up to 100 Km s Does this observation indicate actual mass 
motions or is there another explanation? 
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4. How well do the observed velocity curves represent the 
mass motions of the star? Since the continuous opacity 
scale changes during the pulsation cycle, the observed 
velocities do not necessarily correspond to a given mass 
element. The size of the discrepancy between the observed 
and actual velocity curves is not known. 

5. What is the physical nature of the variable microturbulence 
observed in Cepheids and why is it occasionally supersonic? 

These questions can best be answered by computing hydrodynamic 
model atmospheres. In Chapter II, the method used to compute the 
models is presented. Chapter III contains a discussion of the 
properties of the hydrodynamic envelopes including an investigation 
of the phase lag and second bump. In the next two chapters the 
models are treated as stellar atmospheres; Chapter IV contains a 
discussion of the continuous spectrum; Chapter V, a discussion of 
the line spectrum. The results obtained from these calculations are 
used in Chapter VI to examine several methods used to find Cepheid 
radii. Chapter VII summarizes the results and contains suggestions 


for further work 



CHAPTER II 


METHOD OF COMPUTATION 


A. Diffusion approximation models* 

1. Differential and difference equations . 

The hydrodynamic envelope of a Cepheid can be represented in a 
Lagrangian coordinate system by the following differential equations 


3 dm 


V* 


(II-l) 


¥t + 4 " r2 


d (P + q) * -Gm 
dm 


(II-2) 


d_E 

St 

L - 


= -dL 
a m 

-256 o tt ^ 


(p + q) dJL. 
a t 

r 4 T 3 a_T 
K , a m ? 


(H-3) 

(II-4) 


d r 

Tt " v 


( 11 - 5 ) 


representing conservation of mass (II-l), momentum (II-2), and 
energy (II-3) , energy transport by radiation (II-4) , and the 
definition of velocity (II-5). The symbols used are defined in the 
Notation section. Convection and thermonuclear energy generation 
have been neglected, and the diffusion approximation for radiative 
transfer has been used. The boundary conditions are 


* Most of the material in this section is from working notes prepared 
by G. S. Kutter and W. M. Sparks. Kutter and Sparks (1972) contains a 
summary of these notes. 
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base of envelope: v=0 s L^L^, r=r o (II -6) 

surface: P=P Q , T 4 » | T 4 ff (T+ |) (II-7) 

where Lq and are constants, is the pressure on the surface, and 
T is the Ross eland mean optical depth. When the constitutive equations 
defining P, q, E, and K as functions of T and V are included, the 
differential equations represent a well-defined mathematical 
problem. 

This system of coupled, nonlinear, first-order, partial differen- 
tial equations cannot be solved analytically. To solve the system of 
equations numerically, the star is divided into N concentric mass 
shells whose interfaces are indexed from 1 at the base of the envelope 
to N+l at the surface. These mass shells define a Lagrangian grid, 

^ 1 - / Mq, where - M(l+S ) and 6 is a small number (typically 

-12 

10 ) used to avoid logarithmic singularities at the surface 

(Kippenhahn, Weigert, and Hofmeister 1967). The time grid is defined 
by a set of time steps At 1 *^ = t n+ ^ - t n . In order to decrease 

the amount of interpolation to be done, the variables v, L, and r are 
defined at the interfaces while V, T, P, q, E, and K are defined at 
the midpoints . 

n+l 

The difference equations at time t are 



(II-8) 
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V. 

i 


v n , 
1 - ^ 


. Q 2 , n + % 
+ 6 r i ( q i + % 
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(II-9) 
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+ 0 


(H-12) 


All quantities without a time superscript represent values at time 

n+l 

t 

The parameter 0 allows the expression of the time derivatives 
as a mixture of forward and backward differences. If 9 = 1, the 
equations are implicit; if 6 * 0, the equations are explicit. 

Using 0 > 0 relaxes the restriction put on the time step by the 
Courant condition. 6=1 was used to generate the stable model 
since the time step could become arbitrarily large as the model 
approached equilibrium. Using 0 < 0.5 in equations (II-8) to (11-12) 
produces numerical instabilities as shown by Kutter and Sparks (1972). 
The boundary conditions at the base of the envelope are 

V 1 ' ■ °. WV r l =r 0 (H-13) 

At the surface they are 


N+l 


a - An 


r N+l < P o ~ VlV 


GM 


M o ' Q N+1 ~ V* 


N+l 


(11-14) 
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B N+1 
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r N+l T lH-% 
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^Vkl 
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~ Q n+k ) 
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(11-15) 


The constitutive equations are 

P *= b f T/V + - gas + radiation pressure (11-16) 

3 

E = :rb*T + E + E.. ~ thermal. + radiation (11-17) 

Z r i+e 

+ ionization + excitation energy 

where b* and E., , as well as K are tabulated as functions of T and 
i+e 

V and depend on the composition. In the diffusion approximation 

1 4 

the radiation pressure P^ = j aT , and the radiation energy 
4 

density - aT . 

2. Solving the difference equations . 

The first step in solving the difference equations is to 

linearize by replacing the set jv^, B^, R., j j with the 

setjv 1 +6 V 1 , B. +BB., R. +5R., W i+% +5 W. +% , Z. +% +6 Z. +% ] . 

2 2 

For example r. is replaced by r.(l+26 R ) and P.j, by 

1 i i 



dN 




+ 


dZ 



6Z 


i4% 


The increments are now treated as the unknowns of the system. This 
leads to a system of 5K algebraic equations in 5N + 5 unknowns 
which combined with the 5 boundary conditions, 5 *» 6 * 

fcRj, = 6 W N -f-3/ 2 * s 5Z N+3/2 ~ °* iS a we H" defined system of equations. 
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An initial guess to the model at time t is obtained by 
extrapolating from the preceding model. In general, the difference 
equations will not be satisfied by the extrapolated variables. 

Henyey, Forbes, and Gould (1964) have described a method for solving 
the linearized equations. The system of equations for interface 2 
is set up using the determinant elements given in Appendix A, and a 
Gauss- Jordan reduction is performed to reduce this block to a nearly 
diagonal form. This procedure is continued until the block for inter- 
face N+l has been reduced. A schematic of the partially reduced 
matrix is shown in figure II- 1. Since only the inhomogeneous terms 
and those elements denoted ,1 X n in the figure need be saved, less 
computer memory is required than with other methods. 

After the matrix has been reduced to the form of Figure II-l, 
a back solution is performed to' find the set of increments to be 
added to the first guess. After updating the variables, the above 
procedure is repeated until either the increments or the inhomogeneous 
terms are sufficiently small. 

Although the convergence is nearly quadratic for these models, 
it can be accelerated. If the convergence is monotonic, the 
increments can be multiplied by a number A > 1; if the convergence 
is oscillatory, A < 1 can be used. In these calculations A —0.9 


produced the most rapid convergence. 


When interpolations are needed, geometric means are used, i.e., 
T^ - . T^ ^) . The opacity, though, requires special handling. 

In the hydrogen ionization region (HIR.) the opacity differs by 
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several orders of magnitude between neighboring zones. If the 
geometric mean is used to find the opacity at the interface between 
these two zones, large variations are observed in the light curve every 
time a zone moves through the HIR* According to Stobie (1969a) this 
effect is caused by using too large an opacity at the interface, 
resulting in the zone on the high opacity side of the HIR not radia- 
ting as efficiently as it should. The excess energy retained by this 
zone is released in a very short time when the zone cools and its 
opacity drops. This zoning effect can amount to 0^2 on the light 
curve . 

To avoid this problem, the opacity at the interface is defined 

by 

b l-b 

K ± = K + K* (11-18) 

where K is the larger opacity, K is the lower opacity, and 

T " 

0 < b < 1 is a free parameter. Tests were made to determine the 
best value for b. With b s 0.5, the light curve was very jagged. 

1 

With b < 0.2, the pulsation amplitude decreased, A value of b = - 
did not noticeably affect the pulsation amplitude and produced a 
light curve with bumps due to the zoning of less than O^OS. This 
value of b was used throughout. 

4. Artificial viscosity pressure . 

The artificial viscosity pressure is an arbitrary quantity used 

/ 

to limit the discontinuity at shock fronts (Richtmeyer and Morton 
1967). If no aritifical vi-scosity is included, shock fronts become 

smaller than one zone, and the time steps taken by the model become 

/ 

very short. The artificial viscosity pressure smooths the shock front 
over several zones and thereby limits the discontinuity of the shock front. 
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Kutter and Sparks (1972) have shown that the artificial viscosity 
pressure does not contribute to the mechanical flux and, therefore, 
does not affect the conservation of energy, 

Stobie (1969a) used a constant value of q^ (see equation II- 10) 
throughout the envelope and has shown that the pulsation amplitude 
of his models decreases as the value of q Q increases. Shocks only 
appear in the outer zones, though, and no artificial viscosity is 
needed deeper than the Hell ionization region. A variable q Q can be 
used to account for this fact, i.e.. 


<Vi-% " q e exp ( ' P U /P H>- 

where q is a free parameter of order unity, and P^ was chosen as 
the pressure in the HIR in the equilibrium model. Tests were run 
to determine the best value of q • With q < 2, shocks were large 
and the time step small; with q e >8, the pulsation amplitude decreased. 
A value q =4 was used and resulted in reasonably large time steps 
without affecting the pulsation amplitude. In fact, the pulsation 
amplitude was nearly independent of q fi for 2<q e <6. Changing 
also had little effect on the pulsation as long as q^ was small in 

the Hell ionization region. 

The properties of the models computed using the methods 


presented in this section will be discussed in Chapter III and 
compared to the radiative transfer models discussed below. 

B. Radiative transfer models . 

1. Assumptions . 

The diffusion approximation (equation II-4) represents the 
radiation field very well if the gas Is optically thick. In the 
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atmosphere, however, the radiation field is not isotropic, and a 
more exact solution to the transfer equation must be used. In 
order to simplify the calculations, the following assumptions were 
made: 

1. The atmosphere of the model is plane parallel. Bohm- 
Vitense (1972) has* shown that the effects of spherical 
symmetry are important for yellow supergiants only if 

m bol < ' 7 '* 

2. The gas is in local thermodynamic equilibrium (LTE) . 

This assumption is difficult to justify for a gas as 
tenuous as that of a Cepheid atmosphere but the problem 
would be intractable without it. The results obtained by 
Bell and Rodgers (1967), Parsons (1971b), and Schmidt 
(1971a, b) using LTE indicate that the non-LTE effects are 
not large. 

3. The radiation field can be characterized by one frequency 
group whose opacity is the Rosseland mean. This gray approxi- 
mation has been shown by Bendt and Davis (1971) to differ only 
slightly from the multi- frequency group solution although it 
may lead to an underestimate of the radiation pressure. 

4. Terms of order v/c in the transfer equation can be ignored. 
Since velocities in Cepheids are typically less than 60 km s \ 
these terms should be negligible. 

5. Energy transport is by radiation only. Although Cepheids are 
cool, the density in the envelopes is so low that convection 
is inefficient. In any case, an adequate theory of time- 
dependent convection is not available. 



IS 


With these assumptions, the transfer equation can be written 


as 


p, ~~ « I - ciT 2 * 4 /rr (11-19) 

where p. - cos 9, 
and dT = Kdm. 

The physical parameters needed to solve the hydrodynamic equations 
(XI-1) to (XI-5) are the radiation energy density, E^, the luminosity, 

L, and the radiation pressure, P r , which are given by 



( 11 - 20 ) 

( 11 - 21 ) 

( 11 - 22 ) 


2. Method of solution . 

Due to the anisotropy of the radiation field, the transfer 
equation cannot be directly included in the. hydrodynamic equations 
without introducing a great deal of complexity. It is possible to 
avoid this complexity by noting that the Henyey method is a special 
case of Newton's method for solving systems of nonlinear equations. 
Newton's method can be expressed as 


X (k) = X (k_1) - J-hto), 


(11-23) 


(k) (k-1) 

where x is the vector of unknowns (note that x v - x v is 


the set of increments discussed in Section II-A) , j|(x) is the 
vector of inhomogeneous terms of the linearized equations, and J 
is the Jacobian of _f with respect to x, i.e., “ df^/dx^., and the 

superscript k counts the number of iterations. The system of equations 



19 


is considered solved when | |_f(x)| |< e, where € is the convergence 
criterion. 

If the Jacobian is exact, Newton’s method converges quad- 
ratically (Ortega 1972); if J is only approximate, Newton’s method 
may still converge but the rate of convergence will be slower than 
quadratic. This fact suggests a method for including the transfer 
equation in a Henyey type code. If the code is currently using 
x (k) = x^ k ^ - J^_fp(x) , where D signifies that the diffusion 

(k) (k-1) 

approximation was used to compute the quantity, try x = x s 

j”^f_(x), where R denotes that the transfer equation was used. The 
D 

convergence properties of this method are discussed below, 

3. Convergence of the modified Henyey method . 

Theorem: Given 

(kH) - ..(k) T -U , (k) 


1. x 


= - V4<* w ), k=0,l,2, 

where 

<«K _ 


(11-24) 


2. R = J - J where J is the Jacobian of fL with respect to -x, 
R D R R 


3- 4<x V 0 = 0. 

Show that for all e > 0, their exists a 6 > 0 and an N> 0 
such that for \\ R || < 6, ||x^ - x^| |< e for k > N. In other 

(k) (od) 

words show that x >-*• x as k-» #, 

Proof: 

Equation (11-24) can be expressed as a Picard iteration, i.e., 

x (k+1) = £(x <k) ), (11-25) 

where g(x)= x -. J“^f„(x). By Ostrowski’s theorem (Ortega 1972) 

— — D K 

(co \ 

the series (11-25) converges to x if p (J ) < 1 when J is evaluated 

8 8 

at and if |jx^^ - x^ ^ \\ I s sufficiently small. Here J is the 

6 

Jacobian of £(x) and p (J ) is the spectral radius (i*e. , the largest 

S 

eigenvalue) of J . The j component of £ is 

8 j ■ • 



Therefore J « I - *’I ~ + R ) — R > 

and p (J g ) "P (J“V> * I! ^ R II * II J D X " * U R li‘ 

Choosing 6 = j | J ^ | | ^ implies that p (J ) ^ j | R | j / 6 < 1 if and 
only if | | R| | -1 | Jjj 1 ) ! < I- 

The iteration, therefore, converges if | | R| j < | | J *| | 

Q. E. D. 

This method is analogous to using the secant method as an approxi- 
mation to the Newton-Raphson iteration in the one dimensional case. 
Ortega (1972) has shown that the rate of convergence depends on the 

value of p (J ) . If p(J ) = 0, the convergence is quadratic. If 
S 8 

p ( J ) > 0, as it is when R is non-zero, the convergence will be 
8 

slower than quadratic. Since the matrix R is not known, it is not 
possible to verify a priori whether the iteration will converge to 
the desired solution. Ortega (1972) has shown that, if the procedure 
converges for || x^^ - ^ \\ sufficiently small, it has converged 

to the correct solution. Unfortunately, it is impossible to define 
how small is sufficiently small. The radius of convergence can be 
found only by trial and error. In this code, the first guess is kept 
within the radius of convergence by limiting, the time step. 

4, Tests of the modified Henyey method . 


As a check on this method, a model computed using gray, plane- 
parallel radiative transfer was allowed to evolve to hydrostatic equilibrium. 
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In the optically thick zones, the difference between the temperatures 
in the diffusion approximation and radiative transfer models was about 
0.02%. In the optically thin zones, the model had a T(t) that agreed 
with the exact solution of the gray transfer problem to better than 
0.5%. A shock wave was then sent through the atmosphere. In the 
optically thin zones, a temperature inversion was observed at the 
shock front. No such temperature inversion was observed when the 
diffusion approximation was used. 

Using without modification produced very slow convergence. 

In the diffusion approximation the radiation field at a given point 
depends only on the local state of the gas. When using the trans- 
fer equation, however, the local radiation field depends on the state 
of the gas throughout the model. In a very optically thin zone, the 
radiation field is nearly independent of the local state of the gas. 

To account for this fact,. was modified by multiplying such terms 
as 3B/dT by (1 - e T ). When these modifications were made to J^, 
the rate of convergence improved without affecting the final results. 
The coefficients of the increments used for these models are given 
in Appendix B. 

Using this method with gray, plane parallel radiative transfer 
takes about 30 minutes on the IBM 360/91 to compute one period of a 
Cepheid with a period of 12 days. About 25% of this time is used 
to solve the transfer equation. It would, therefore, take 8 additional 
minutes per frequency point for the non-gray problem, roughly a factor 
of four faster than the method of Keller and Mutschlecner (1970) and 
comparable in speed to the variable Eddington approximation used by 
Bendt and Davis (1971). 
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This method can be easily modified to include such effects 
as spherical symmetry or non-gray radiative transfer by replacing 
the transfer subroutine used here. In fact, the method has wider 
applicability. For example, a Henyey type evolution code could be 
modified to include deviations from hydrostatic equilibrium or 
convective overshoot. 

This modified Henyey method has been used to compute deep 
envelope models of a 12 d classical Cepheid. The following chapters 
contain a discussion of the properties of these models. 
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Figure XI- 1. Schematic of the partially reduced matrix of the 

coefficients. The x’s denote non^zero matrix elements. 
All elements not shown are zero. 


chapter III 


ENVELOPE OF THE MODEL 


A. Equilibrium model . 

1. Selection of the model . 

There are two important criteria used to select a model for a 
study such as this. First, the model should represent a well- 
observed star to allow a comparison of theory with observation. 

Since there are bright Cepheids of many periods, this criterion does 
not greatly restrict the selection of model parameters. Second, the 
model should not present unnecessary computational problems. Stobie 
(1969c) has shown that short period Cepheids are more likely to be 
pulsating in a harmonic than in the fundamental mode and some may be 
pulsating in two modes simultaneously. Long period Cepheids also 
present problems. These stars are cool enough that convection is 
expected to be important. In addition, there is evidence of strong 
shocks in their atmospheres (Rodgers and Bell 1968; Dawe 1969). 

These two criteria indicate that the model should be chosen 
to have a period between 8 d and 15 d . One further criterion is the 
availability of other theoretical results. Both Keller and 
Mutschlecner (1970, 1971) and Bendt and Davis (1971) computed models 
with periods of about 12 d . For these reasons a model with a period 
of -12 d was selected having T e ^ f =5700°K, L^SOOOI^, and M=5F^(Stobie 
1969c). 
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2. The First Guess 

Once the model has been selected, the first guess to the model 
at t = 0 must be generated. Although experience has shown that the 
radius of convergence of the Henyey method is large, a reasonably 
accurate first guess to the initial model is needed. Simply taking 
a published envelope model will not suffice since difference in 
zoning, opacity, and mean molecular weight will cause problems. 

Either the iteration for the first time step will not converge or 
the model will develop such strong shocks that it will eject a large 
fracticnof its mass. To avoid these problems the static envelope code 
used by Rose and Smith (1970, 1972) was modified to include the King 
IVa tables for opacity and mean molecular weight (King 1972) . Table 
IIX-1 contains the number fractions of the elements included in cal- 
culating these tables. 

After a static envelope having the desired mass, luminosity, and 
effective temperature was generated, the zoning criteria were selected. 
The total mass of the envelope was chosen so that the base of the en- 
velope had T^IO 6 °K since Stobie (1969a) has shown that including 
deeper zones does not appreciably affect the pulsation properties of 

the model. In order to have enough optically thin zones to study the 

26 

atmosphere, a mass for the first zone of 1.5 x 10 gm was selected. 

It is also necessary to have enough zones in the model to perform the 
numerical integration over depth with sufficient accuracy. The 
results of Stobie (196ya) and Bendt and Davis (1971) indicate that at 
least 50 zones are required. In order to minimize any problems associ- 
ated with the zoning, 100 zones were used. This combination of 100 

26 

of the surface zone of 1.5 x 10 gm, and mass of the 


zones, mass 
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envelope of 0.5 resulted in a ratio of the mass of neighboring 
zones of 1,15. Once these parameters were selected, it was possible 
to interpolate in the static envelope model to obtain a first guess 
to the model at t=0. 

3, Relaxation to equilibrium . 

The first guess to the model at t-0 was used as input to the 
hydrodynamic code. Since the static envelope code included 
convection and the hydrodynamic code did not, the model was not 
quite in hydrostatic or thermal equilibrium. In order to decrease 
the computer time required for the model to evolve to equilibrium, 
the pulsation was artificially damped by converting some of the kinetic 
energy into thermal energy at each time step. With this damping the 
model reached equilibrium in about 4000 time steps. 

The equilibrium model was found to have a slightly different 
effective temperature than the static envelope model. Some adjustment 
was, therefore, necessary. The two free parameters that determine the 
position of the model in the H-R diagram are the luminosity incident on 
the lower boundary, L^, and the radius of the lower boundary, r^. In 
equilibrium the luminosity leaving the surface must be equal to 
Thus, if Lq is the luminosity selected for the static model, the 
effective temperature of the model can be changed by varying r^. If 
is decreased, the model becomes more compressed, and increases; 

if r^ is increased, T decreases. After changing r^ the model must 
be forced to return to equilibrium. This procedure of changing 
and forcing the model to return to equilibrium was continued until the 
model had the desired effective temperature. 
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The final equilibrium model, STB, has L-5000 1^, M-5^, T ^ * 

5730°K, and R^-71.7 Rq. STB has 25 optically thin zones ranging from 
T = 4.6 x 10 ^ to t = 3. Due to the large opacity in the hydrogen ioniza- 
tion region, the first optically thick zone is at T = 580. It is 
possible, therefore, to refer to the optically thin zones as the atmo- 
sphere without ambiguity. Since the atmosphere of STB has a thickness 
of only 1 % of the radius of the model, the plane parallel assumption 
used in the radiative transfer calculations is justified. 

Table I II -2 contains the parameters of the final equilibrium 
envelope, and figures III-l to I1I-3 show the logarithm of the tem- 
perature, pressure, and density, respectively, versus the logarithm of 
Lagrangian mass coordinate. Note the large temperature increase and 
the density inversion in the hydrogen ionization region. Both of 
these effects would be smaller, but not eliminated, if convection 
were included in the model (Latour 1970; Eoll 1973). However, the time 
it takes a convective element to travel one mixing length is comparable 
to the period of the Cepheid. Convection cannot be included correctly 
in the hydrodynamic models until a theory of time dependent convection 
is developed. 

B. Full amplitude model . 

1, Growth to full amplitude . 

Once a static envelope with the desired luminosity and effective 
temperature has been generated, the pulsation can be initiated in one 
of two ways, by soft self-excitation or by hard self -excitation 
(Ledoux and Walraven 1958). The former approach allows the pulsation 
to grow from the computer round-off "noise 11 (Cox ana Giuli 1968, 
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pg 1125). The soft self-excitation requires no prior knowledge of 
the pulsation mode or the full amplitude velocity distribution, but 
the model must be followed for 5 to 10 e- folding times. The 
second approach has been used by Christy (1968) and Stobie (1969a) and 
consists of imposing a finite amplitude velocity distribution on the 
static model. The initial velocity field is usually chosen so that 
the initial kinetic energy amplitude is about 25 % of its full amplitude 
value. If the imposed velocity distribution does not correspond to 
nearly a pure mode, transients will have a large amplitude and will 
take several e- folding times to die out. This method assumes a sub- 
stantial amount of information about the full amplitude pulsation, but 
it allows full amplitude to be reached in only 2 or 3 e-folding times. 
Since Stobie (1969a) gives the initial velocity distribution for a wide 
variety of Cepheid models, and because e-folding times of Cepheids are 
of the order of 100 periods, the hard self -excitation approach was 
used. 


Since the radiative transfer models take about five times more 
computing time per period than the diffusion approximation models, 
the diffusion approximation was used in generating the full amplitude 
model. The initial velocity distribution was given by 


where 

and 


v = -v 0 (r / r p ) 6 , 

r •= radius at T « 1, 
P 

Vq “ 10 km s 


The maximum of the total kinetic energy in the envelope during a period, 
KE max » provides the most reliable measure of the growth of the pulsation. 
After initiating the pulsation, it was found that KE 

max 


decreased for 



29 


tWo periods as the high order harmonics introduced by the initial 

velocity distribution died out. After about 8 periods these transients 

had a much lower amplitude than the fundamental mode. 

In order to speed the growth to full amplitude, an artificial 

amplification was used (Stobie 1969a). Choosing a phase when all 

zones were moving outward with nearly their maximum velocity, the 

velocity in each zone was multiplied by 1.5, doubling the kinetic 

energy. Although this procedure introduced some transients, they 

quickly died out. At period 26 the method of opacity averaging at 

the HIR was changed to the method described in Chapter XI. This 

change greatly reduced the zoning effects resulting in fewer shocks 

and, therefore, less damping due to the artificial viscosity. The 

rate of increase of KE then increased. The model was then allowed 

max 

to pulsate for another 120 periods, about 3 e- folding times (Stobie 

1969c). At this point KE was increasing by only 0.02% per period. 

' max 

Period 150 was chosen as the full amplitude, diffusion approximation 

model, IHC (for implicit hydrodynamic code) . Figure III-4 shows the 

approach of to full amplitude. 

IHC was used as the starting model for the radiative transfer 

calculation. As can be seen in figure 111-4, KE max decreased smoothly 

(indicating few transients) and approached a value about 10% lower 

than KE in IHC. After 30 more periods KE „ was decreasing by less 
max max 

/ 

than 0.1% per period. Period 180 was chosen as the full amplitude 
radiative transfer model, RET. The smaller amplitude of RDT relative 
to IHC is due to changes in the structure of, and, therefore, the 
work done by, zones passing through the HIR. The ratio of /PdV in 
RDT to JPd V in IHC for those zones passing through the HIR is 0.85 
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accounting for the lower amplitude of RDT. 

2. Properties of the full amplitude models . 

The period of RDT is 12?05, which combined with the mas a and 
radius of STB gives Q » p/p /pQ 0.046. Christy (1968) has shown 
that Q=G. 022 (R /Rq) / M)^. Using the radius and mass of STB 

gives Q = 0.043 in good agreement with RDT. 

A comparison between the full amplitude RDT and IHC models 
shows that, aside from differences due to the lower amplitude of 
RDT, the two models produce nearly identical light and velocity 
curves conf inning the results of Bendt and Davis (1971) and Davis 
(1971). The major difference between the models is that IHC never 
has temperature inversions in the atmosphere while RDT does whenever 
there are sufficiently strong shocks present. 

It is not difficult to understand this difference between the 
models. Consider a stellar atmosphere in radiative equilibrium. 

In the diffusion approximation 


L ** dB — ®i+l , 

dT T i+ i * T i 

If the temperature in zone i+1 is increased, the luminosity will 
increase by an amount 


AB 


Ah 


i+1 


T . , - - T. 

1+1 1 


If T ^ is small, as it is in the atmosphere, AL will be large. 

The excess thermal energy in zone i+1 will be radiated away and the 

atmosphere will rapidly return to equilibrium. On the other hand, 

the luminosity given by the transfer equation is 

T . 


I ~ J B(t)E 2 (t-T)dt - f 1 B(t)E 2 (T i -t)dt, 
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where E 2 is the second exponential integral. Perturbing the temperature 


in zone i+1 results in 

AL ~ Ab 


i+1 


r 


E„(t-T.)dt. 
2' x 


If t - T . is small, 
x+1 i 




< T i+ rV- 


In the atmosphere, then, AL will be small and the return to equilibrium 
will be slow. This behavior of the radiative transfer solution has 
been discussed by Whitney (1967). 

Since the temperature inversions in RBT occur only above t*= 0.01, 
they will be important only for spectral lines and will have only a 
small effect on the continuum. In fact, the continuum forming regions 
of the atmosphere are nearly in radiative equilibrium at nearly all 
phases as predicted by Whitney (1967) . In the following discussion all 
values will be taken from RDT, but the conclusions reached apply 
equally well to IHC. 

The relative radius change 

AR _ R MAX ' _ n , , 

R “ 'ha** ' ‘ 


is in good agreement with observed values (Nikolov and Tsvetkov 

1972), as is the velocity amplitude v -v . - 45 km sec ^ (Stibbs 

max mm 

1955). Figure III- 5 shows the bolometric light curve. The small 
ripples on the light curve are due to the zoning effects discussed in 
the preceding chapter and are small enough not to confuse the inter- 
pretation of the main features. The asymmetry of the light curve 
falls within the observed range (Nikolov 1968). In contrast, the 
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models of Keller and Mutchlecner (1970, 1971) and Bendt and Davis 
<1971) are too asymmetric in the sense that the rising branches of 
their light curves are too steep. Since they use a shallow envelope 
driven from below by a coarsely zoned deep envelope model, and RDT 
was computed using fine zoning throughout, this difference is 
probably due to zoning. 

Aside from the overall variation of the light curve and the 
zoning effects, there are three features of interest, a sudden 
decrease in light output at phase $ = 0.15, a shoulder on the rising 
branch near $ *= 0.2, and a bump on the falling branch near $ = 0.6, 

All three features show zoning effects superimposed on them and, 
therefore, are probably not due to the zoning. The shoulder on the 
rising branch and the bump on the falling branch will be discussed 
in the next section but it should be noted that Cepheids with periods 
near 10^ often show similar features (van Genderen 1970). 

The dip at $ - 0.15 is disturbing since it is not observed 
(Nikolov 1968) . Bendt and Davis (1971) also see a similar feature 
on their light curves as do King, Cox, Eilers, and Davey (1973) in 
their coarsely zoned diffusion approximation model of an 11^5 Cepheid. 
Since this feature is usually dismissed as a zoning effect, an 
attemp.t was made to find its cause. A detailed search of the computer 
output revealed nothing related to the zoning that could cause the dip. 
Hillendahl (1968) attributes the feature to the artificial viscosity. 
There is a pressure wave moving outward that produces the shoulder 
at $ & 0.2. The artificial viscosity causes the pressure to rise 
ahead of the temperature. The pressure rise causes the density to in- 
crease which, in turn, raises the opacity. When this region of in- 
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creased opacity reaches t~ 1, the amount of light emitted decreases. 

An opacity increase of 3% is sufficient to produce the dip observed 
on the theoretical light curve. The same mechanism will produce 
a dip preceding the bump at ? 3 0.6. The bump on the falling branch 
would look more like a shoulder if the artificial viscosity could be 
removed. 

Figure III-6 shows the velocity corresponding to T = 0.2. The 

shoulder on the rising branch cannot be seen, but the bump on the 

falling branch is quite pronounced. Since the velocity curves are 

usually used to classify the bumps in the models, this model would be 

described as having a single bump on the falling branch. According 

to the Hertzsprung sequence a Cepheid with a period of 12 should have 

a bump on the rising branch. If the velocity curve was used to 

classify this model, it would be concluded that the model had the 

wrong mass. This problem of classifying bumps on the light curves 

can contribute to the Cepheid mass discrepancy (see Chapter I). 

-3 -1 

Figure III-7 shows the velocity curves for T - 10 ,10 ,1 

and for a mass zone having T “ 0.2 in STB* The feature on the 
rising branch appears only for the most optically thin zones and 
may be responsible for the peculiarities in the cores of the Ca II 
H and K lines observed at this phase (Kraft 1967). The progressive 
nature of the wave travelling through the atmosphere is readily 
apparent. Substantial velocity gradients are present in the atmo- 
sphere from § - 0.3 to*! = 0.7. The velocity curve following a given 
mass element is the same as that for T - 1 except in the vicinity of 
the second bump even though t a i moves through more than 10 mass 
zones. This point will be considered further in Chapter VI. 
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3, Phase lag . 

The phase lag between the light and velocity curves, A $, 
depends on which velocity curve is used. The phase shift ranges 
from A $ = 0.06 for TaltoA^ B 0 for T * 10 . Linear, 

adiabatic theory predicts that light maximum should occur at minimum 
radius, not near equilibrium radius as observed. Explanations for 
the phase lag range from Eddington's suggestion in 1917 that the 
phase lag is a natural consequence of the processes limiting the 
amplitude of the star to Christy's in 1968 which attributes the phase 
lag to a skewing of the light curve due to non-linear effects. Since 
the linear calculations of King, Cox, Eilers, and Davey (1973) show 
the phase lag, it is probable that the lag is the result of non- 
adiabatic, not non-linear effects. 

Castor (1968) treating the HIR as a discontinuity, has sug- 
gested that a theory proposed by Eddington (1926) is correct. 

The large heat capacity of the HIR delays light maximum. Since 
the HIR lies at the top of the transition region between the quasi- 
adiabatic envelope and the non-adiabatic atmosphere, the luminosity 
gets "frozen- in" at the top of the HIR. 

The phase lag can be seen in Figure III-8 which is a 3-D plot 
showing the variation of luminosity, L/Lg^, as a function of mass 
point and phase as viewed from the center of the star. Note that 
phase increases from right to left. The inset is a schematic repre- 
sentation and will be used to define points of reference in the figure. 
Point A is in the He II ionization zone. This part of the model is 


nearly adiabatic, and, as expected, the luminosity maximum coincides 
with radius minimum. By the time the He I ionization zone is reached 
at point B, there is a substantial phase shift. A further, small phase 
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shift is introduced in the HIR, the region between B and C. The 
"freezing- in” of the flux in the atmosphere, point C, is seen as 
luminosity perturbations moving outward at constant phase. 

Plate I is a different representation of the same data and shows 
two periods of the motion. The abscissa is the time coordinate, the 
ordinate is the Lagrangian radial coordinate. Large values of L/L Stb 
appear as bright areas while small values appear dark. The base of 
the envelope is at the top of the figure. The bright area nearest 
the top of the figure corresponds to point A in figure III-8. As can 
be seen the phase of light maximum increases continuously until the 
atmosphere is reached. At this point, the flux becomes "frozen- in" and 
the maximum moves outward at constant phase. The model, therefore, 
suggests that the HIR plays only a small role in generating the phase 
Shift. The phase shift appears to vary continuously through the 
transition region between the quasiadiabatic envelope and the non- 
adiabatic atmosphere. 

Plate I shows another interesting phenomenon. In the He II ioniza- 
tion region light minimum follows light maximum by about half the period. 
In the He I ionization zone this phase difference is still nearly half 
the period. Only in the atmosphere is the light curve very asymmetric. 
The asymmetry of the light curve is, therefore, either an atmospheric 
phenomenon or due to the HIR. This effect will be discussed in the 
next section. 

C.- Cause of the second bumps . 

1. Christy mechanism . 

In 1926, Hertzsprung classified Cepheids by the shape of their 
light curves (Payne-Gaposhkin 1951). Short period Cepheids have very 
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asymmetric, smooth light curves. As the period increases the curves 
become more symmetric and a second bump appears on the falling branch 
in Cepheids with a period near 7^. This bump moves closer to the 
primary maximum until at 10^ the bump appears to coincide with maximum 
light. Near a period of 12^ the light curves often show two shoulders, 
one on the rising and one on the falling branch. Cepheids with longer 
periods show a single bump on the rising branch. Near a 15^ period 
the curves again become smooth and asymmetric. Figure III-9 taken 
from Payne-Gaposhkin (1951) illustrates this sequence. 

Christy (1970) using the results of his nonlinear calculations 
explains the second bump in the following manner. Near the phase of 
maximum compression the He II ionization zone is rapidly heating 
causing the zone to expand. This expansion sends a pressure wave out- 
ward which appears at light maximum. The expansion also sends a pres- 
sure wave inward which is reflected by the core and appears at the 
surface during the next cycle as a second bump. He then finds that 
the time delay, D (in days), can be used to find the radius of the 
model from R / Rq = 4.05D. Since the period of the star depends on 
both its mass and radius, Christy can find both the mass and radius 
of a star from its period and the phase of its bump. In general, 
masses found in this way are about half those predicted by stellar 
evolution theory. 

It should be noted, however, that the Hertzsprung sequence is an 
average property of Cepheids. Van Genderen (1970) has shown that the 
phase of the second bump of individual Cepheids having the same period 
varies over a wide range. He also shows ttet, while the Hertzsprung 
sequence holds on the average up to a period of 10^, there is almost 
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no correlation between bump phase and period from 10 to 30 . 

Beyond 30^ the bump appears near minimum light. In those Cepheids 
having two bumps on their light curves, he suggests that two inde- 
pendent phenomena are operating. 

2, Hillendahl's mechanism . 

The mechanism proposed by Christy explains part of the Hertz- 
sprung sequence. It does not, however, explain the double bump 
Cepheids. Another mechanism is needed. Figure III- 10 shows the 
velocity as a function of mass point and phase® The inset will be 
used to define points of interest in the figure. Effects due to the 
zoning have been labelled "z” and are quite small compared to the 
main features. The shaded area in the inset shows the bottom side 
of the surface while the dashed line follows the HIR. The point 
corresponding to maximum light is labelled f? A n ; B is the shoulder on 
the rising branch of the light curve, and C is the bump, on the falling 
branch. The line marked D is an inward moving pressure wave discussed 
below, and line E indicates the location of the inward moving pressure 
wave described by Christy (1970). 

Figure III-ll is a different view of the same data. After reaching 
maximum expansion velocity near D the atmosphere begins to slow down 
under the influence of gravity along line E. However, a disturbance 
originating at point A changes the sign of the acceleration and propo- 
gates both outward toward B and inward toward C. The velocity reaches 
a second maximum near F and then decreases under the influence of 
gravity along G. The curve marked H indicates the velocity curve 
deeper in the envelope. Figure III-12 shows the origin of the inward 
and outward moving pressure waves. 


The line marked A is the locus of 
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points of maximum velocity in the atmosphere; line B marks the locus 
of points of maximum velocity in the envelope. As is evident, the 

atmosphere reaches maximum velocity before the envelope and starts 
to slow down while the envelope is still accelerating compressing a 
region near the HIR. This compression generates pressure waves moving 
both outward (line AB on figure III-ll) and inward (line AC on figure 
III-ll) . Line D(E) shows the inward (outward) moving pressure dis- 
turbance described by Christy (1970). 

Figure III- 13 clearly shows that the Christy pressure wave from 
the preceding cycle (line C) arrives at the surface of the model near 
light maximum (A), not near the second bump (B) . In fact, this 
pressure wave reaches t * 1 very near = 0.2 and is responsible for 
the shoulder on the rising branch of the light curve. 

Figure III-14 illustrates one possible explanation for the second 
bump proposed by Hillendahl (1969, 1970), The points A, B, and C 
correspond to the three local maxima on the light curve at phases $ = 0.2, 
0.35, 0.6, respectively. The Christy pressure wave reaches the top of 
the envelope at point D and is responsible for the feature at B, 

According to Hillendahl, the local velocity minimum at E is the result 
of a rarefaction wave moving inward. The feature at C can then be 
attributed to a secondary, "blow- off" shock. The inward moving 
pressure disturbance (line C in figure III-12) is then a second rare- 
faction wave. In Hillendahl 's interpretation, the features labelled 
z in figure III-10 are further blow-off shocks. This effect can also 
be seen in plate II. Here large positive velocities appear bright, 
zero velocity gray, and large negative velocities dark. Two periods 
are shown, i.e., phase increases from 0 to 2 from left to right, and 
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the base of the envelope is at the top of the figure. The local 
velocity minimum preceding the second bump (line AB in figure III-ll) 
can be followed backward in time into the envelope. Hillendahl 
considers this minimum to be a result of the rarefaction wave follow- 
ing the deep envelope pressure disturbance (line C in figure IIX-13) . 

The inward moving pressure wave (line AC in figure III-ll) he associates 
with a second rarefaction. 

There are several problems with this interpretation. Hillendahl 
predicts as many as 5 blow-off shocks per period. All features on 
the light and velocity curves following the bump on the falling branch 
can be associated with zones moving through the HIR. While these 
zoning effects may be masking the secondary blow-off shocks, it is 
unlikely that the shocks would not be seen at all. Another difficulty 
is the inward moving pressure wave (AC in figure III-ll). If this 
feature is to be associated with rarefaction following the secondary 
shock causing the bump, it should follow the local velocity maximum. 

It does not. It appears to originate before the second velocity maxi- 
mum. The Christy mechanism cannot explain this feature either. Hillen- 
dahl's mechanism also predicts that all Cepheids should have a second 
bump on the falling branch since the primary expansion should always 
cause a secondary, blow-off shock. 

3. Atmospheric oscillation mechanism . 

Since neither Christy's nor Hillendahl's mechanisms adequately ex- 
plain the Hertzsprung sequence a third mechanism was sought. Inspection 
of figure III- II shows that the surface layers appear to be pulsating 
nearly sinusoidally from point J to point B with a period roughly 2/3 
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the period of the envelope. Figure HI-15 shows the power spectra 
of two zones obtained from Fourier transforms of their velocity 
curves. The dashed line refers to a zone at the top of the envelope. 

No secondary periodicities are apparent with less than 3 times the 
frequency associated with the envelope period, The amplitude of 

the feature near V = ifO is so low that it may be due to noise genera- 
ted in the transform process or to the inward moving pressure wave. 

The solid line, referring to an atmospheric zone, indicates a secondary 
periodicity near with about 1/3 the amplitude of the main pulsation. 
The higher frequency features are probably the result of zoning effects 
or noise generated in the transform process. 

The secondary periodicity can be understood in the following 
analysis which was used by Christy (1962) in an early attempt to 
predict the properties of Cepheid atmospheres. According to Lamb 
(1932) the critical period of an isothermal atmosphere is given by 

F a = v g /(yg), (XXX- 1) 


where 



y ~ 5/3 for a neutral gas, 

_gm_ M 


Substituting values taken from STB shows that 

d (R/Ry ) 2 

P = 0.0025 - 7 - 7 --^; 

a (M/l^) 


which gives = 2?5 for STB. This period is too small to account 
for the secondary periodicity seen in the atmosphere of the model. 



41 


Gough, Ostriker, and Stobie (1965), on the other hand, have 
shown that a better approximation to the Cepheid atmosphere is an 
atmosphere with a constant temperature lapse rate. Choosing values 
from STB results in 


d (R ' V 
p a = °-oo7 omg 


(III-3) 


or P =7° which is about 0.6 times the period of RDT. It appears 
a 

that the bump on the light curve of RDT occurs because the atmospheric 
period, P a , is comparable to, but less than, the envelope period, P^. 

The occurrence of two bumps on Cepheid light curves can be ex- 
plained qualitatively as follows: Cepheids with periods less than 7 d 

have P < 0.5 P . Since the atmosphere is being driven far from 
a e 

resonance, the amplitude of the atmospheric mode is low, and the 

atmosphere follows the motion of the envelope. Starting at 7 d , the 

driving frequency approaches the resonant frequency of the atmosphere, 

but not until about 10 d is the amplitude of the atmospheric mode large 

enough to produce an observable shock from the compression of the HIR. 

The multiple bumps seen on the rising branch of ultra-violet light 

curves of £ Doradus by Hutchinson (1974) may indicate that the atmo- 

sphere is beginning to produce these shocks. In the period range 7 

to 10 the bumps appearing on the falling branch are produced by the 

Christy mechanism. From 10 d to 12 d the Christy bump appears on the 

rising branch and the atmospheric oscillation bump appears on the 

falling branch. As thp period increases beyond 12 d the amplitude of 

the atmospheric mode grows, but the compression of the HIR decreases 

as the atmosphere and envelope begin to oscillate in phase. There 

d d 

are no bumps from 15 to 30 since the compression of the HIR is too 



42 


small to produce observable shocks. At about 25 = P g , the 

atmosphere is oscillating at its maximum amplitude, but there are no 
bumps since the atmosphere and envelope are always in phase. Beyond 
30^ bumps appear near minimum light as the envelope begins to slow 
down before the atmosphere leading to compression of the HIR near 
minimum light. 

As shown by Payne -Gaposhkin (1951) Cepheids with periods near 
10 d are anomalous compared to those at 8 and 12 . There is a dip in 
the velocity amplitude-period relation and in the upper envelope of 
the light amplitude-period relation near P = 10^ . The anomaly can be • 

understood if primarily envelope oscillations are being observed up 

d d 

to 8 , and atmospheric oscillations are being observed beyond 12 . 

4. Ratio of the atmosphere to the envelope period . 

Further calculations will be needed to check the validity of the 
above picture, but it can be examined for consistency using published 
relations among L, P , R, and M. In the following discussion, R and 
M are given in solar units and periods in days. In logarithmic form, 
equation (III-3) becomes 

log P * -2.14 + 2 log R - log M (III-4) 

while the relation V ^ ]p Q = Q(M,R) results in 

log P « a + b log R + c log M. (III-5) 

6 

Combining the period -luminosity law and the mass- luminosity relation 
results in an equation of the form 
log M = d log + e. 


(III-6) 
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where the constants d and e depend on the choice of a mass -luminosity 
relation. The scatter introduced into (II1-6) by multiple crossings 
of the instability strip will be small since 90% of the observed 
Cepheids are expected to be in the second crossing (Iben 1966). 
Similarly, Fernie (1968) has shown that a relationship of the form 
log R - f log P e + g (III-7) 

can be expected. If equations (III-6) and (III-7) are to be 
consistent with (II1-5) 

a + bg + ce = 0, (III-8) 

bf + cd “ 1 

The values of the constants in equations (III-5) to (III-7) 
can be found in the literature. Christy (1970) and Fricke, 

Stobie, and Strittmatter (1972) give essentially the same values 
for the constants in (III-5), namely 

log » -1.62 4- 1.72 log R - 0.68 log M. (I1I-9) 

Using Stobie 's (1969c) mass -luminosity law gives 

log M = 0.34 log P e + 0.33, . (III-10) 

where the constant, e, has been adjusted slightly to give the correct 
mass for RDT . Fernie (1968) has used the Wesselink method to find 
log R = 0.56 log + 1.24, (XIX-11) 

where the zero point has been adjusted downward as suggested by 
Parsons (1972). 

Equations (III-9) to (III-ll) do not satisfy the consistency 
conditions (III-8) . If the period-radius law, which was derived 

directly from observations, is redefined so the consistency conditions 
are satisfied 

log R ** 0.7-2 log P + 1.07. 

e 


(III-12) 
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The large change in the coefficients in the period-radius relation 
indicates that either the wrong mass-luminosity law was used, or 
there are systematic errors in the radii determined using the 
Wesselink method. The former possibility will be investigated below, 
and the latter in Chapter VI. Equation (III-12) agrees almost ex- 
actly with the relationship given by Wooley and Carter (1973) indi- 
cating an error in the Wesselink method as used by Fernie and Parsons. 


The adopted relationships are 

log P - -1.62 4- 1.72 log R - 0.68 log M, (III-13) 

log M « 0.34 log P e + 0.33, (III-14) 

log R * 0.72 log P e + 1.07 (III-15) 

log P a ® -2.14 + 2 log R - log M. (111-16) 

Substituting equations (111-14) and (III-15) into (HI-16) gives 

log P = 1.10 log P - 0.33, 
a e 

N 

or log (P a /P e ) - 0.10 log P e - 0.33 (111-17) 


Equation (III- 17) indicates P = P near log P =3 which does not 

a e e 

support the explanation of the double bump Cepheids give above. 

Repeating this analysis using the results of stellar evolution 
theory as given by Iben and Tuggle (1972a) and again defining the 
period- radius law from the consistency conditions gives 

log P £ - -1.53 + 1.73 log R - 0.79 log M, (III-18) 

log M = 0.30 log F + 0.56, (111-19) 

log R *= 0.72 log P^ *f 1.14. (III-20) 

Note that the period- adius law (III-20) is nearly the same as 
(III-15). The discrepancy with the observed relation (III-ll) is, 
therefore, not due to the choice of a particular mass- luminosity law. 
Substituting (III-19) and (III-20)- into (III-16) gives 
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log (P B /P ) = 0.21 log P e - 0.30, (111-21) 

indicating P s P near P = 30 . 

One of two conclusions can be reached from this analysis. Either 

the analysis of the atmospheric oscillations is in error, and the 

bump on the falling branch of the light curve has a different origin, 

or the mass-luminosity relation predicted by evolution theory is 

more nearly correct than that of pulsation theory. There is other 

evidence to support the evolution masses. If masses determined from 

the phase of the second bump are ignored, both Fricke, Stobie, and 

Strittmatter (1971, 1972) and Iben and Tuggle (1972, a,b) can explain 

the mass discrepancy. Since the pulsation masses depend critically 

on the calibration of the observations, reasonably small changes in 

the T vs. (B-V) relation, the helium abundance, and the zero point 
ef f 

of the period-luminosity law can remove the mass discrepancy. But, 
if there are two mechanisms which can produce bumps on .the light curves, 
the bump masses should be ignored. 

Although this discussion favors the evolution masses, the 
question cannot be settled without computing more models with many 
optically thin zones, a time consuming procedure. Additionally, 
more accurate values for the constants in equations (111-5) to (111-7) 
are needed. In particular, the variation of the pulsation constant, 

Q, with mass and radius is not well-known. The results of Cogan (1970), 
Cox, King, and Stellinwerf (1972) are considerably different from 
those of Christy (1970), Fricke, Stobie, and Strittmatter (1972), and 
Parsons and Bouw (1971). Further study including the effects of 


convection is needed. 



TABLE III-l 


NUMBER ABUNDANCES USED 
FOR KING IVa TABLES 

X = 0.70 Y * 0.28 


ELEMENT 


abundance 


H 

He 

C 

N 

0 

Ne 

Na 

Mg 

A1 

Si 

Ar 

Fe 


9.07156E-1 
9.13793E-2 
2.84443E-4 
8.01673E-5 
6 .36793E-4 
3 .58809E-4 
1.42560E-6 
1.79473E-5 
1. 18576E-6 
2 .27738E-5 
2 .38399E-5 
3.69354E-5 
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Figure III-l. log T VS, log (1-Q) in the equilibrium model 
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Figure III-3. log/Jvs. log (1-Q) in the equilibrium model 
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Figure 111-6 ♦ Velocity at t = 0.2 vs. phase for the full amplitude model 
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Figure I II -8. Luminosity vs. mass point and phase. Every second mass 

point has been plotted. Point A is in the Hell ionization 
zone; point B , the top of the quasi-adiabatic envelope-, 
point C, the atmosphere. 
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Figure Ill -9. Observed light curves of Cephcids illustrating the Hertzsprung 
sequence, (from Payne-Gaposhkin 1951). 
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Figure III-14. Velocity vs. mass point and phase 
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FREQUENCY (1/PERIOD) 

Figure 111-15. Power spectra of velocity curves. atmospheric 

zone , envelope zone . 



CHAPTER IV 


CONTINUOUS SPECTRUM 

A. Atmosphere calculation . 

1, Snapshot approach . 

In the preceding chapters the model has been discussed in terms 
of pulsation theory. Since many optically thin zones were included 
in the calculation of the hydrodynamic models, it is also possible to 
discuss. the model in terms of stellar atmospheres. The hydrodynamic 
equations (II-l) to (II-5) are most easily solved if a Lagrangian 
mass coordinate is used as the independent variable. On the other 
hand, the light emerging from a model is more easily computed if the 
optical depth, T, is the independent variable. The physical parameters 
of primary interest are also different. For example, while the total 
pressure is needed to calculate the hydrodynamic motions, the electron 
pressure is more important for determining the emergent spectrum. 

To facilitate the stellar atmosphere calculations, the optically 
thin zones of the hydrodynamic models were used as snapshots of the 
atmospheric structure. Since all relevant physical variables change 
on a time scale much greater than the time it takes a photon to diffuse 
through the atmosphere, no large errors are introduced by using this 
snapshot approximation. This assumption is the same one made in 
dropping terms of order v/c in the transfer equation. After describing 
the methods used to convert the hydrodynamic models to a form suitable 
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for the atmosphere calculations, the properties of the continuous 
spectrum will be discussed in the remainder of this chapter. Chapter V 
will be devoted to a discussion of the line spectrum. 

2, Converting the models . 

The hydrodynamic models at 200 phases at equal phase intervals 
of 0.005 period were converted to a form with Rosseland mean optical 
depth, dr = Kpdr, as the independent variable. The atmospheric models 
were computed starting at the surface and were continued inward until 
the temperature reached 25,000°K. In all cases the base of the 
atmosphere was at a large enough optical depth that including more 
points would not change the computed spectrum. 

Before computing the emergent spectrum, the electron pressure, 

P^, had to be determined. Since P^ depends on the number of free 
electrons determined from the Saha equation, and is needed to 
compute the ionization fractions, an iteration must be. performed to 
find P^. This iteration is described in Appendix C. Knowing the 
variation of temperature, density, radius, and electron pressure with 
optical depth makes it possible to compute the spectral energy dis- 
tributions of the models. The monochromatic opacity, K^was computed 
with a code written by Bell (1974) that includes the continuous 
opacities of H, H , He, He , Sil, Mgl and Rayleigh scattering. 

Next the monochromatic optical depth scale was computed from 



and the transfer equation was solved to give the flux at the desired 
wavelength, . 



The wavelength region from 912JL to Sp. was divided into 21 
wavelength bands whose boundaries occur at the absorption edges of 
the continuous opacity sources, A Gauss-Lobatto quadrature using 
4, 6, or 10 wavelength points was performed for each band. The 
sura of the contributions of these bands represents the total flux 
of the model 


-/ 


oo 134 

V* " i w i ■ CT W" 


(IV- 1) 


Equation (IV- 1) was used to define the effective temperature of the 
model, T ef ^ s at each phase. 

The total flux was defined as the integral of not F since 

a V 


,-2 


10 


-1 


F^(\ * 5p.) / max 10 while F v (X - 5p) / max 

Since there is essentially no flux shortward of 912j[ for stars in 
this temperature range, the integration can be computed with fewer 

/ oo 

B^dA with these 
u 

quadrature points resulted in an error of 0.3%. Before these fluxes 
can be used to compute colors, the effect of the Bpectral lines must 
be taken into account, 

3. Line blocking approximation . 

In the range of spectral types populated by Cepheids, the contribu- 
tion of spectral lines to the total opacity is large in certain wave- 
length regions. Because including the effect of lines in the mono- 
chromatic opacity would require an excessive amount of computer time, 
the line blocking approximation was used. In this approximation, the 
flux in each wavelength band is multiplied by a line blocking factor, 

0 £ T[ £ 1, which represents the fraction of the flux absorbed by the 


lines in the band. The line blocked flux, F„, is then defined by 

B 
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where F y is the flux computed including only continuum opacity sources. 
This approximation does not conserve energy since the flux removed 
by the lines is simply ignored; but, since gray radiative transfer was 
used to compute the hydrodynamic structure, it was felt that use of 
the line blocking approximation was warranted. 

Figure IV-1 shows TJvs.X for a model taken from Parsons (1969) 
having = 6000°K and log g » 1.8. The line blocking factors are 

those of Bell (1974) . Figure IV-2 illustrates the difference between 
F„ and F T1 . The dashed line is the spectral energy distribution 
computed using only continuous opacity sources; the solid line in- 
cludes the line blocking factors shown in Figure XV- 1. Note the 
linear scale of the ordinate chosen to accentuate the difference 

between f and F„. It is clear from Figure IV-2 that the line blocking 
B U 

has a reasonably large effect on the U, B, and V magnitudes but changes 

the R and I magnitudes only slightly. 

The line blocking factors are in the form of tables of T) vs, X 

for a set of T -- and g Before the line blocking factors can be 

eff eff 

applied to the hydrodynamic models, T ^ and g e ££ must be defined. 

The method used is described in the following section. 

4. Effective temperature and gravity . 

The parameters of major importance in the hydrodynamic calculations 
were the luminosity emerging from the model and the radius of the photo- 
sphere. The line blocking factors, though, are tabulated in terms of 

T and g ■ The effective temperature of the model can be found by 
ext err 

2 4 . 

solving L »» 4 ttR CTT e £f» but determining is more complicated. A 

stellar atmosphere is hydrostatic equilibrium has a unique gravity. 
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g = GM / r . The hydrodynamic models, on the other hand, are 
experiencing accelerations and each zone has an effective gravity 
given by 


8 


ef f 


- dP =- GM + r, 

dm 2 

r 


(IV-3) 


where r is the acceleration of the zone. A suitable method for 
finding the mean value of l n the line forming region of the 

atmosphere was needed. 

A stellar atmosphere is normally defined in terms of and 

& e ff ^ e ff * s a measure of the temperature of the model; of 

the pressure. For this reason, the mean effective gravity of the 

Cl . 

atmosphere, & e ff> was defined from equation (IV-3) as 
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r* I dP 

_ JO \ dm 


eff 


-t/to . 
e dT 


£ 


V t/to at 


(IV-4) 


where - 0.1, was chosen to 'limit the average to the -line forming 
region. This value of g ef £ - was used along with to determine 

the line blocking factors at each phase. 

The variations of T e ^ and with phase are shown in Figures IV-3 

and IV-4, respectively. The large spike in log S eff and the dip in 
T g ££ near phase § =0.15 may be artifacts of the artificial viscosity 
(see Chapter III). The large increase in log g e £ f near § = 0.5, on 
the other hand, is real and is caused by a shock moving out through 
the atmosphere. If the feature' near $ = 0.15 on the log g e ££ curve is 
ignored, the variation is found to be less than has been observed by 
Schwartzchild, Schwartzchild, and Adams (1948), Schmidt (1971b), and 
Parsons (1971 a,b). This difference is not significant since g eff of 
RDT was found by averaging over the entire atmosphere while the observed 



dynamic gravities refer to a small part of the atmosphere. The g ef £ 
variation of a given zone is much larger but spurious accelerations 
due to zoning effects are also enhanced. 

Oke, Giver, and Searle (1962) in their analysis of the RR Lyra 
star SU Draconis found a range of phases having roughly constant 
^eff 3n< * ®eff’ As s ^ own * n Fi 6 ure IV-5, RJDT behaves differently. 
Ignoring the effects of the atmospheric snock wave and the spike 
caused by the artificial viscosity (by following the dashed lines in 
Figure IV-5) gives an open curve in the log S e ££ ~ log T e ^£ diagram. 
During rising and falling light, the atmosphere moves at nearly 
constant S e ££> i«e., isobarically, while near maximum and minimum 
light, it changes isothermally. If it is assumed that the pressure 
at some point in the atmosphere is proportional to g & ££» and that the 
perfect gas law with T g ^£ and this pressure can be used to find the 
specific volume, V, a P - V diagram can be constructed (Figure IV-6) . 
The dashed lines in Figure IV-5 have been followed in constructing 
this diagram. Following the actual gravity variation instead results 
in the addition of two long, narrow loops containing little enclosed 
area. Since these loops needlessly clutter the diagram and contribute 
little to the discussion, they are not shown. The area enclosed by 
this curve indicates that the atmosphere has a destablizing effect on 

the model. Although this curve is not a proper thermodynamic integral 

/ 

since it does not follow a given mass element, Christy (1962), who 
constructed a similar diagram from observations of SU Dra, concludes 
that this P - V diagram reflects the influence of the HIR on the star. 
5. Defining the color system . 

Once T’ -- and g have been found for each model, the line 
eff erx 

blocking factors can be applied to the spectral energy distributions 



computed above, and the broad band UBVRI colors can then be calculated. 
The relative filter sensitivity functions for the U, B, and V filters 
were taken from Azusienis and Straizys (1969) and those of R and I from 
Johnson (1964). These magnitudes were on an arbitrary scale. In order 
to make comparisons with observations, V, U~B, B-V, V-R, and R-I were 
converted to Johnson's system using least squares fits to 

Cj - a + bC, (IV-5) 

where 

C = color in Johnson' system, 

C - color in arbitrary system. 

Line blocked spectral energy distributions of Parsons' (1969) models 
were used for the calibration. Unfortunately, there was no single 
source available that gives the Johnson colors needed for the trans- 
formations (IV-5). The (U-B) and (B-V) T colors were taken from Bell 

and Parsons (1974); (R-I) was taken from Schmidt (1973); (V-R) _ was 

J J 

computed from the calibration by Caputo and Natta (1973). The absolute 
visual magnitudes, V, were computed from the bolometric magnitude using 
the bolometric correction of Kraft (1961) . All colors given in the 
following are on Johnson's system. 

The absolute visual magnitude light curve, vs. $ is shown in 
Figure IV-7. Although the zoning effects are smaller than 0^05, it 
was found that they masked features of interest on the HR and color- 
color diagrams. The light curves were, therefore, smoothed to minim- 
ize these zoning effects. These smoothed light curves shown in Figure 
IV-8 were used for the comparison with observations discussed in the 


next section. 



B , Comparison with observations . 


1. Light curve parameters . 

The light curves computed in the previous section were compared 
with observations. Since the model was not selected to represent a 
specific star, the theoretical results were compared with the mean 
properties of Cepheids with periods near 12 d . In Figure IV-8, the 
increase of light amplitude with decreasing effective wavelength is 
readily apparent. Somewhat less apparent, but still discernable, 
is the phase shift noted by Stebbins, Kron, and Smith (1952) and 
Wisniewski and Johnson (1968). This phase shift is in the sense that 
long wavelength light curves lag short wavelength light curves. The 
phase shift between the U and I curves of RDT is about 0.04 period in 
good agreement with observations. Table IV*- 1 compares the light ampli- 
tude of RDT with observed values from Wisniewski and Johnson (1968) . 
Unfortunately, the only Cepheid they observed with a period near 12 
is £ Gem, a star with a low-amplitude, symmetric light curve. The 
table shows that the RDT amplitudes fall within the observed ranges. 

RDT can also be compared with the light curve parameters published 
by Schaltenbrand and Tammann (1971). Table IV-2 compares RDT with all 
classical Cepheids in the period range 11^ to 13^ that they observed. 
Again the values for RDT fall within the range of observed quantities. 
Figure IV-9 and Table IV-3 give the adopted light and color curves 
for RDT. 

2. Color-temperature relations . 

A further check on the models can be made. By treating each of 
the 200 models as an independent observation, the color- temperature 
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laws can be computed. Figure IV- 10 shows the dependence of (U-B), 
(B-V), (V-R), and (R-I) on log T e ^. The dependence of (U-B) on 
g e ^ produces the open curve in the (U-B)-log plane. The 

(B-V), (V-R), and (R-I) curves, on the other hand, closely approxi- 
mate straight lines and should, therefore, be good temperature 
indicators. The fits to the relationships 

log T eff - a + b(X-Y), 


6 eff “ c + d(X-Y), 

where (X-Y) represents the color, are given in Table IV-3, as are 
the relationships derived by Schmidt (1971a). For comparison, 

Kraft (1961) gives log T g ££ - 3.886 - 0.175 (B-V) and Rodgers (1970) 
gives 0 = 0.64 + 0.337 (B-V). 

This agreement with the observed values is much better than 
expected and results from partial cancellation of two errors. 

Ignoring convection produces a model having too large a temperature 
gradient in the photosphere. Using the line blocking approximation, 
which neglects the backwarming effect of the line opacity, produces 
a model having too small a temperature gradient. Ignoring both 
effects results in a temperature structure closely approximating that 
of the stars studied. Since the observed col or -temperature laws were 
derived using nonpulsating super-giants and hydrostatic model atmospheres, 
the accuracy of the color temperature fits indicates that deviations 
from radiative equilibrium must be small, as predicted by Whitney 
(1967). The use of static model atmospheres to study the colors of 
Cepheids should, therefore, produce reliable results. 



C, Loops in the (U-B)-(B-V) diagram . 

Having established that the colors of RDT reproduce the observed 
colors allows the models to be used to answer questions requiring a 
theoretical approach. One of these questions has already been 
answered. The hydrodynamic atmospheres are nearly in radiative 
equilibrium at most phases and can be approximated by a series of 
hydrostatic model atmospheres. 

Another interesting question is the cause of the loops in the 
(U-B) - (B-V) diagram, Abt (1959) has suggested that the loops are 
caused by 

a) excess ultraviolet emission from shocks, 

b) sensitivity of the continuous opacity, K , 
to the electron pressure, P , 

c) unusual line blanketing, 

d) lines being partially filled in by ; emission lines, 

e) continuous emission possibly originating in a chromosphere, or 

f) non-thermal dependence of line strength (i.e., with P ) . 

Figure IV- 11a shows the color- color diagram for RDT. The curve 
resembles that of T) Aql and would be classified by Uikolov and 
Kunchev (1969) as being linear or nearly linear. The curve is notice- 
ably open from phase $ = 0.1 to § = 0.5 having a maximum width of 0*1*04 
in (U-B). Since neither emission lines nor chromospheric emission was 
included in the model, the loops could not be produced by (d) or (e) . 
Figure IV- lib shows the color -color diagram for RDT excluding the line 
blocking factors. Due to the change of scale the curve appears to be 
much more open but still has a maximum width of 0.05 in (U-B). Since 
lines have not been included in calculating this case, the openness 



75 


cannot be explained by (c) or (f). There is a strong shock in the 

atmosphere only near $ = oTs and the excess emission amounts to about 
m 

0.02 in (U-B) . Possibility (a), therefore, can be excluded. Only (b) , 

the sensitivity of K to P » remains. In the temperature range of RDT, 

v> e 

the continuous opacities of H and H are the primary opacity sources. 

In the atmosphere, most of the free electrons come from the metals, 
aid the ionization of the metals depends linearly on. P^ through the 
Saha equation. Since the wavelength dependences of the H and H 
continuous opacities differ, a change of P^ at fixed temperature will 
change H /H and, therefore, the wavelength dependence of K , Since 
the U filter contains the Balmer jump, this effect will be more pro- 
nounced in U than in B or V, producing a loop in the color-color diagram. 

D. Mean colors of Cepheids . 

Figure IV-12 indicates another problem. The model traces an 
open path in the HR Diagram running nearly parallel to the lines of 
constant period and covering the entire width of the instability strip. 
Many attempts at defining period- luminosity or period- luminosity-color 
relations depend on the mean values of the luminosity and colors. 

There are at least 3 distinct methods for computing these means, 

a) the intensity mean of the color-(B-V) , 

b) the magnitude mean of the color-(B-V) , 

c) the difference of the intensity means of two 

magnitudes -{ b)^. - (v) , 


Unfortunately, these method do not give the same results. Table 
IV-4 compares these means with the colors of STB. Method (c) is 
normally considered to be the most physically meaningful average 



since conservation of energy in a steady state system requires 



where L is the bolometric luminosity and P the period. There is 
nothing to require this relationship to hold in a limited wavelength 
region, however, since the wavelength dependence of the luminosity 
varies during the cycle. Table IV-4 indicates, however, that method 
(c) does best reproduce the values for STB confirming the results of 
Cox and Wing (1973). The discrepancy between the mean values of 
(U-B) from methods (a) and (b) and that of STB is difficult to 
explain, but there appear to be two possibilities. First, since the 
U filter contains the Balmer jump, it is sensitive to the confluence 
of the hydrogen lines near 3650ft. and, therefore, to g e ^. Second, 
the U filter contains the wavelengths with both the largest and the 
smallest opacities. The U magnitude, therefore, is affected by a 
larger region of the atmosphere than the other filters. Deviations 
from radiative equilibrium will affect the U magnitude for a larger 
fraction of the cycle than it will the other magnitudes. In either case 
the larger amplitude in U will accentuate the differences among the 
averaging schemes. 

E. Zero point of the P-L relation . 

One further point of interest is the location of STB in the HR 
Diagram. The + in Figure IV- 12 marks the location of STB and the 
dashed lines show the Sandage and Taramann (1969) instability strip 
and two lines of constant period. The location of STB indicates 
that the zero point of the Sandage and Tarnmann P-L-C relation should 
be lowered by 0.2. Then and Tucele f"I971b^ suepest that this laree a 

w ' / ' ' — - — >= - 
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change is reasonable considering the uncertainties in defining the 
zero point. Unfortunately, they suggest an upward revision to remove 
the discrepancy between the pulsational and evolutionary masses. On 
the other hand, recent studies of the zero point using secular paral- 
laxes indicate the adopted zero point might be too high. Jung (1970) 
finds that the shift should be AM = 0?4 ± 0^4 while Wielen (1974) 
gives AM = 0^2 ± 0?4. Wielen points out that the increase in the zero 
point required by Iben and Tuggle cannot be ruled out and is just barely 
tolerable. 

Due to the approximations made in computing the models, the results 
presented in this chapter must be considered tentative. More accurate 
models are needed to verify the correctness of these conclusions. The 
inaccuracies of the models are even more important when the line pro- 
files are considered in the next chapter. 
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Table IV- 1 


Comparison of observed and 
computed light amplitudes 


STAR 

P 

AU 

AB 

AV 

Ar 

Al 

X CYG 

16 *4 

2.4 

1.6 

1.0 

0.60 

0.42 

RDT 

12.1 

1.9 

1.6 

1.0 

0.55 

0.45 

U AQL 

7.0 

1.5 

1.2 

0.8 

0,52 

0.42 





Table IV -2 





Comparison of observed and 




computed 

light curve 

parameters 


star 

AV 

A(B-V) 

A(U-B) 

A$ 

RX 

aur 

.692 

.356 

.461 

.505 

SS 

cma 

.992 

.457 

.485 

.553 

XX 

CAR 

.835 

.457 

- 

.642 

RY 

CAS 

.969 

.494 

.528 

.599 

XX 

CEN 

.886 

.453 

.508 

.511 

KK 

CEN 

.979 

.539 

- 

.528 

SU 

CRU 

.678 

.663 

.668 

.664 

FZ 

CYG 

.778 

.525 

.533 

.519 

AA 

GEM 

.672 

.434 

.457 

.532 

UU 

MUS 

1.031 

.582 

.705 

.636 

U 

NOR 

.967 

.472 

.553 

.541 

SY 

NOR 

.890 

.334 

.319 

.646 

SV 

PER 

.787 

.404 

.293 

.609 

Z 

SCT 

1.001 

.549 

.831 

.596 

TY 

SCT 

.881 

.453 

.678 

.534 

DR 

VEL 

.570 

.342 

- 

.544 

RDT 

.99 

.49 

.46 

.57 
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Table IV -3 

Adopted light and color curves 


$ 

B01» 

N 

U-B 

B-V 

V-R 

R-I 

0.00 

“4,051 

-3.922 

0.648 

0.807 

0.674 

0,461 

0.05 

-4.132 

-4.039 

0.553 

0.736 

0.627 

0.427 

0.10 

-4.275 

-4.223 

0.445 

0.628 

0.555 

0.380 

0.15 

-4.347 

-4.310 

0.381 

0.568 

0.517 

0.351 

0.20 

-4.510 

-4.500 

0.345 

0.490 

0.454 

0.314 

0.25 

-4.537 

-4.525 

0.351 

0.506 

0.463 

0.313 

0.30 

6.35 

-4.657 

-4.647 

0.327 

0.439 

0.426 

0.298 

-4.758 

-4.758 

0.338 

0.423 

0.402 

0.289 

0.40 

-4.720 

-4.717 

0.360 

0.460 

0.426 

0.306 

0.45 

-4.612 

-4.592 

0.383 

0.513 

0.470 

0.336 

0.50 

-4.457 

-4.403 

0.434 

0.611 

0.544 

0.372 

0.55 

-4.405 

-4.337 

0.467 

0.642 

0.565 

0.390 

0.60 

-4.421 

-4.354 

0.483 

0.649 

0.567 

0.393 

0.65 

-4.369 

-4.290 

0.543 

0.704 

0.597 

0.411 

0.70 

-4.280 

-4.174 

0.607 

0.763 

0.641 

0.436 

0.75 

-4.195 

-4.062 

0.678 

0.833 

0.681 

0.457 

0.80 

-4.102 

-3.944 

0.717 

0.843 

0.706 

0.482 

0.85 

-4.030 

-3.854 

0.754 

0.871 

0.725 

0.497 

0.90 

-3.972 

-3.782 

0.781 

0.914 

0.746 

0.503 

0.95 

-3.973 

-3.802 

0.752 

0.886 

0.725 

0.493 



Table IV -4 


Color-temperature relations 


108 T eff = 3 + b(X ‘ Y) 

®ef£ = C + d(X ' Y) 

Present Present Schmidt (1971) 



a 

b 

c 

d 

c 

d 

B-V 

3.877 

-.174 

.642 

.349 

.641 

.309 

V-R 

3.906 

-.252 

.583 

.505 

.543 

.593 

R-I 

3.918 

-.394 

.561 

.791 

.574 

.773 


Table IV- 5 

Magnitude and intensity mean colors 



U-B 

B-V 

V-R 

R-I 

V 

<X-Y) I 

.507 

.654 

.571 

.393 


< x ‘ y >m 

.518 

.665 

.576 

.396 

-4.261 

(X> I '<Y> I 

.451 

.614 

.552 

.386 

-4.302 

STB 

.433 

.634 

.575 

.386 

-4.273 
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Figure IV-1. 



Line blocking factors, 
r ~ 6 000 °K ami log g = 


vs . 1 for a model with 

1 n 

A . O . 
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Figure IV-2. Spectral energy distribution for a model with = 6000°K 

and log g = 1.8. no line blocking; 

with line blocking. Note linear scale of ordinate. 



F^(10 7 erg s 







Figure IV-4. log T vs. phase. Mote decrease in at same phase 

as increase in g near $ = 0.15 




LOG geff 



LOG T e ff 

Figure IV-5* log 8 ef £ vs, log T ^ Dashed lines show the curve 
following dashed line in Figure IV-3. 
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Figure IV-6. P~V diagram constructed using T e ^ and The dashed 

lines in Figure IV-3 has been followed in constructing 
this curve. The .numbers indicate phase. 
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Figure IV-7. Absolute visual magnitude, vs. $ 
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Figure IV-8. Smoothed light curves in U, B, V, R, and I 
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Figure IV-9. Adopted light and color curves for full amplitude model 
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Figure XV- II. a) (U-li) vs. (B-V) including line blocking factors 
b) (U-B) vs. (B-V) without line blocking factors 
Numbers indicate phase 
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Figure IV- 12, Curve traced by the model in the H-R diagram. Dashed 
lines indicate the Sandage and Tammann instability 
strip and two lines of constant period. Numbers 
indicate phase. 





CHAPTER V 


LINE SPECTRUM 

A, Method . 

1. Basic approach . 

In the preceding chapter it was shown that the broad band colors 
computed from the hydrodynamic model atmospheres agree with observa- 
tions of Cepheids, Ideally, one would also like to show that the 
models reproduce the observed line spectrum. There are several 
reasons why a detailed comparison was not made. The primary reason 
is the large amount of computer time required to perforin a synthetic 
spectrum analysis. In addition, the models do not include line 
blanketing, non-LTE effects, or non-gray radiative transfer, effects 
that are expected to be more important in the line-forming region 
than in the photosphere (Bohm-Vitense 1972; Osmer 1972). Since the 
number of variables that must be changed to investigate the range 
of observed phenomena is large, a somewhat heuristic approach was 
adopted. A single line, the X4494.57A line of Fel (excitation potential 
2.2eV), was chosen to give a reasonable variation of the number of 
absorbers with depth in the atmosphere. The equivalent oscillator 
strength, f, was treated as a free parameter. By varying f it was 
possible to study the effect of the moving atmosphere on lines of 
different strengths, 

2, Computational procedure* 

The method for computing line profiles in a moving atmosphere 
follows a suggestion made by Chandrasekhar (1945). The line 

* , 

The material in this section is from Karp (1973). 
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absorption is described by a Voigt profile, H(a,u), where a is the 
damping parameter and u = (\-Xq)/AX^ 0 When there are velocities in 
the atmosphere, u must be modified to account for the motions of the 
gas. This can be done by letting 

u' * u + p,vX/c 

where p,v is the local velocity of the gas projected onto the line 
of sight s and using 11(8,^) to compute the line opacity. The 
specific intensity can be computed from 

f co 

^(t^expC-^Ai) dt^/p. , 

and then the flux from 

y»l 

F v (0) =2 / ^(0, p>dp 

*0 

Other methods commonly used, such as the Feautrier method or the 
quadrature integration of Milne's second equation (Kourganoff 1952), 
cannot be used since the opacity is a function of p. 

B, Microturbulence . 

1, Velocity gradient mechanism * 

Before considering the detailed line shapes, the equivalent widths 
of the lines, W, were used to investigate microturbulence in Cepheids. 
Struve (1932) introduced microturbulence to explain the anomalously 
large Doppler broadening velocities he found in supergiants. This 
was explained by Struve and Elvey (1934) as being due to either "a 
turbulence of small eddies 1 ' or to "several shells which expand with 
different velocities." The concept of microturbulence as a turbulence 
of small eddies has been attacked by Worrall and Wilson (1972), They 
claim that the concept of microturbulence is valid only in terms of 

^ ■ r j- s- _ 

The material in this section has been adapted from Karp (1973). 



the reversing layer treatment of line formation and attribute the 
large Doppler broadening velocities to inadequacies in the theory 
of line formation, particularly the LTE assumption. It has also been 
suggested by Evans and Schroeder (1972) and Andersen (1973) that 
systematic errors in the measured equivalent oscillator strengths 
are responsible for microturbulence. Since neither of the above 
suggestions has yet been shown to be responsible for the observed 
microturbulence, there is still a great deal of interest in small 
scale turbulence. Hearn (1974) has recently shown that the flux 
of energy required to maintain microturbulence is at least 100 times 
greater than the acoustic flux generated by convection in the sun. 

In supergiants the problem is not as bad but the acoustic energy is 
still a factor of 10 too small to maintain the observed microturbulence 
even though the acoustic flux is much higher (de Loore 1970) . 

It has been found that micro turbulence varies with height in the 
atmosphere of supergiants (Wright 1946; Huang and Struve 1960) and 
with phase in Cepheids (van Paradijs 1972). Differential motions have 
been reported in the atmospheres of supergiants of many spectral types 
(Abt 1957; Rosendahl and Wegner 1970; Aydin 1972) and, in particular, 
in Cepheids (van Hoof and Deurinck 1952 ; Sanford 1956; Dawe 1969). 

Van Paradijs has noted that the microturbulence in Cepheids is a 
maximum near the phase of most rapid contraction. Dawe (1969) and van 
Hoof and Deurinck (1952) have shown that the velocity gradient is 
appreciable at this phase. Sanford (1956) has shown that line widths 
increase just before maximum light in T Mon and SV Vul, a phase at 
which he observes a large velocity difference between weak and strong 
lines. The effect of a velocity gradient on the curve of growth has 
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been investigated by Abhyankar (1964 a,b), Kubiowski and Ciurla (1965), 
Ciurla (1966) and Arakelyan (1969). Although they showed that a 
velocity gradient raises the flat part of the curve of growth mimicing 
raicroturbulence, no attempt was made to compare computed and observed 
velocity differences. Before investigating the velocity gradient 
mechanism using the hydrodynamic models, a preliminary study was made. 

2. Test of the velocity gradient mechanism . 

To make the test case as realistic as possible, a model atmo- 
sphere with T ff * 6300°K and log g=1.8 from Parsons (1969) and the 
Fel line discussed above were used. A wide range of effects was 
studied by computing curves of growth for log a = -1, -2, -3 and for 
micro turbulent velocities 5=0 and 5 km sec * by varying the number of 
absorbers in the line of sight. These curves are shown in Figures V-l 
to V~3 . The curves for 5=0 and 5=5 km sec * do not come together at 
large because the ordinate is -log W/X instead of -log 

Underhill (1947) has shown that a velocity of expansion (or contraction) 
constant in T cannot change the equivalent width, W, of a line. Such 
a velocity field will produce asymmetric profiles, however, due to the 
integration over the surface. As a check, several profiles were 
computed with v(T) ■ 20 and 40 km s” 1 . In no case was the change in 
W greater than 1 per cent. This change is due to errors in the angle 
and frequency integration and can be used as a rough estimate of the 
numerical errors in all these calculations. 

For these preliminary calculations, an arbitrary choice of v(T) 
was made, v(t) = - a log T. This is convenient because it allows a 
correlation of mean optical depth of formation of a line and its 
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observed radial velocity and is nearly linear with geometric height 
in the line forming region. An arbitrary constant may be added to 
v(T) without changing W, but it will change the shape of the profile. 

The results for §=0 and 5 km s ^ and a=*5 and 10 km s ^ are 
shown in Figures V-l to V-3 . With a c 10 fen s \ the curve of growth 
(C^iq) is nearly identical to the normal curve of growth with 5 km s * 
(Cjjj.) until the damping portion is reached. In all cases C^q has a 
wider plateau than C^. An observer would interpret this as being 
due to a lower value of the damping parameter, a. The decrease in a 
at phases when § is large has been observed by Rodgers and Bell (1964, 
1968). This change is a is easily understood. The vertical shift 
between the damping parts of C^q and is proportional to the ratio 
of the Doppler widths. Since C ^ and C have the same Doppler width 
as C^q, the three curves must join in the strong line asymptotic limit. 
The only way this can happen is for tp be below as the lines 

get strong. 

Dawe (1969) has plotted observed velocity versus mean optical 
depth of formation for weak lines in & Car. Reading from Dawe's 
Figure 3, the weakest lines which are formed near T = 0.3 have a 
velocity of about 18 km s ^ while those formed near T ~ 0.1 have a 
velocity of 22 km s Rodgers and Bell (1968) have observed l Car 
and find that 5 = 7.5 km s ^ at this phase. If the above results can 
be extended to this case, a * 15 km s ^ should be used to correspond 
to 'the observed microturbulence which gives v(0.1) - v(0.3) = 7.2 km s 


Correcting for the integration over the surface by the factor p - v 


v ^ — 1.375 (Parsons 1971) gives a predicted velocity difference of 
5.2 km s" 1 which, considering the uncertainties involved, is in excellent 
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agreement with the 4 km s' 1 velocity difference observed by Dawe. 

The radial velocities observed from the minima of the computed 
profiles are given in Table V-l. An exact comparison is not meaningful 
since the observed velocities and shapes of the lines are more sensitive 
to the velocity distribution than is the shape of the curve of growth. 

3 „ Microturbulence in the hydrodynami c models. 

The hydrodynamic Cepheid atmospheres make it possible to study 
the velocity gradient mechanism using more realistic velocities than 
used in the preceding section. The microturbulence was determined 


for each of 20 models at equal phase intervals, A$ - 0.05, by perform- 
ing a differential curve of growth analysis relative to STB. Figure V-4 
shows the curves of growth for STB computed by varying log gf for the 
selected Fel line. After computing similar curves for the 20 hydro- 
dynamic models, the curves were shifted horizontally until their linear 


parts coincided. The remaining vertical shift is 

XJL 


joBvJ « 


2 2 
1 + v 


log 


STB 


where DBV stands for the Doppler broadening velocity, v^ and Vg^g 
are kinetic velocities of Fe atoms in the hydrodynamic and stable 
models, respectively, and l is the unknown microturbulent velocity. 

If it is assumed that the kinetic temperatures of the two models are 
in the ratio of their effective temperatures, 5/v gTB can be found. 

The kinetic velocity in STB was found from [rav] measured in Figure V-4 

ef 

In Table V-2, which summarizes the results, Alog W - DBV , is 
the ratio of the Doppler broadening velocities, v^/v gTB is the ratio 
of the kinetic velocities, and 5 is the microturbulent velocity, and 



AV is discussed below. All velocities are in km s \ 

It is clear that the velocity gradients are not responsible 
for all the microturbulence observed in Cepheids. Typical micro- 
turbulent velocities in Cepheids are 5 to 10 km s ^ (Rodgers and 
Bell 1964, 1968; Schmidt 1971b; van Paradijs 1972). The latter 
value is supersonic. They also find that § varies with phase with 
an amplitude of 2 to 5 km s ^ . The microturbulence in RDT only has 
an amplitude of 1 km s ^ . Thus, the velocity gradient mechanism 
cannot explain the variable microturbulence unless more accurate 
models show that the effect has been underestimated here. However, 
since | rarely exceeds the sound speed by more than 1 or 2 km s ^ 


in Cepheids, the velocity gradient may be responsible for the super- 
sonic part of the microturbulence. 

Since it appears that two phenomena contribute to the observed 
microturbulence, an attempt was made to see how the classical micro- 
turbulence, 5 and the microturbulence caused by the velocity 

gradients, § , combine to produce the total observed microturbulence, 

8 r 

A test case was computed using the model at $ - 0 including 

£ , - 2 km s ^ in addition to § - 1.4 km s It was found that 

cl gr 

= 3,1 km s 1 indicating the Doppler broadening velocity is consistent 
with 


DBV = /Hi + (t +? ) 2 . 

V cl gr' 

Since the total microturbulence is not the square root of the sum of 
the squares but is more likely even a relatively small 


§ can make § appear to be supersonic, 
gr c 



Schmidt (1974) has attempted to correlate the observed micro- 
turbulence with the observed velocity difference between weak and 
strong lines and concludes that there is no apparent correlation. 

Even when the velocity distribution in the atmosphere is known, 
as it is in the models, specifying the correct velocity difference, 

A V, for such a correlation is difficult. Since the velocity 
gradient is not montonic, it is not clear how an average A V charac- 
teristic of the atmosphere can be found. Figure V-5 shows § versus 
A V where A V was defined as the maximum minus the minimum velocity 
on 'the interval 0.01 < t < 1. The correlation is weak and there is 

a good deal of scatter, making it unlikely that the correlation 

© 

could be found observationally . The line shown, a least squares fit 
to the data, is given by 

c = 0.79 + 0.21 AV , 

The line does not pass through the origin due to the problem defining 
A V mentioned above. 

From this analysis, it appears that microturbulence in Cepheids 
is still a mystery. The velocity gradient mechanism is not significant 
in producing the observed micro turbulence and may not be large enough 
to explain the observed variation of the microturbulence with phase. 
Further work is needed. 

C. Line profiles . 

1. Asymmetries . 

Having computed hydrodynamic model atmospheres makes it possible 
to compare theoretical and observed line shapes. There are two 
mechanisms which contribute to the asymmetry of the line profiles. 



One of these is a geometric distortion first studied by Chandrasekhar 
(1945). The observed flux profile is a weighted sum of intensity pro- 
files from different parts of the stellar disk. As illustrated in 
Figure V-6, each intensity profile is Doppler- shifted by an amount, 
p,Vp, where p, - cos 0 and - pulsation velocity of the atmosphere. 

One such set of intensity profiles is shown in Figure V-7 and illus- 
trates an important computational point. Enough angle points must 
be included in the angle quadrature to insure sufficient overlap of 
the intensity profiles. If too few points are used, spurious bumps 
appear on the flux profile. Expressing the total width of the in- 
tensity profile in velocity units, v^, gives a condition on the maxi- 
mum angle separation, Ap. = v^/v. Sharper lines and higher velocities 
require more angle points. 

Even a constant velocity in the atmosphere will produce asymmetric 
line profiles. According to van Hoof and Deurinck (1952) a line will 
appear symmetric unless 


V_ > 2.5 (AX) 

c H 

where is the rest wavelength of the line; v, the velocity in the 
atmosphere; and (AX)^, the half width of the undisturbed line. Thus, 
weak lines will appear more asymmetric than strong lines. Figure V-8 
illustrates this point. The weak line, (a), is very asymmetric while 
the strong line, (b), appears undistorted even though both lines were 
computed with a constant atmospheric velocity of 20 km s \ The 
dashed line in (a) is the bisector of the weak line and is characteris- 
tic of lines distorted by this geometric effect. 



105 


The second mechanism distorting the line profile is most noticeable 
in strong lines. Bell and Rodgers (1964) and Kraft (1967) have suggested 
that the asymmetry of these strong lines is due to velocity gradients 
in the atmosphere, i.e., different parts of the line are formed in 
layers of gas moving with different velocities. Figure V-9 illustrates 
the velocity gradient asymmetry. Not only is this line too strong to 
show a noticeable geometric distortion, but its bisector has a shape 
very different from that in Figure V-8(a). The profile in Figure V-9 
can be compared to the 45081 line of Fell in P Dor observed by Bell 
and Rodgers (1964). 

The geometric distortion is expected to be most important at 
phases of maximum velocity while the distortion of strong lines should 
be greatest at phases of maximum velocity gradient. In view of the 
difficulty of determining the velocity gradient from the velocity 
difference between weak and strong lines, it appears the asymmetry of 
strong lines is a better indicator of the velocity gradient in the 
atmosphere. 

2. Cheshire cat lines . 

Another interesting phenomenon is the occurrence of "Cheshire 
■Cat" (Carroll 1865) lines. This phrase was used by Underhill (see 
Kraft 1967, pg. 240) to describe the extra component of strong lines 
often observed in Cepheids. These extra components seem to have no 
antecedents, but suddenly appear at phases when strong shocks are 
expected in the atmosphere. 

Grenfell and Wallerstein (1969) have attempted to explain the 

splitting of the line in SV Vul in the following way. The red 

• -I 

component of H is due to gas falling inward at 70 km s , 

CG 


which 
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appears to move a distance of half a stellar radius. The violet com- 
ponent, which has a velocity 40 km s ^ more negative than the metal 
lines, is associated with the pulsation of the atmosphere. When 
these two layers of gas collide, emission lines of H and He should 
be observed but are not. Wallerstein (1972) suggests that the inward 
moving material originates in a circumstellar shell. Skalafuris 
(1974) suggests a different mechanism based on the work of Whitney 
(1956), but also asserts that the splitting is due to velocity differ- 
ences. 

Kraft (1967) has stated that the velocity difference between the 
two components may be as large as 30 to 40 km s ^ for the low exci- 
tation metal lines while Grenfell and Wallerstein (1969) and Wallerstein 
(1972) report velocity differences of up to 100 km s ^ in the H lines 

of SV Vul and T Mon. Since the atmosphere of RDT never has velocities 

-1 -1 
greater than 30 km s or velocity differences greater than 15 km s , 

the occurrence of the "Cheshire Cat" line shown in Figure V-10 indi- 
cates a different origin for the secondary component. The deeper com- 
ponent has a velocity of 19 km s ^ while the shallower has *3 km s \ 
a difference of 22 km s Examining the model shows the maximum 
velocity difference is 3 km s * and the mean velocity is 10 km s \ 

The splitting of the line is obviously not caused by differential motions. 

Further inspection of the model revealed that there was a tempera- 
ture inversion of about 300°K near T - 3 x 10 , the region in which 
the line core is formed. To verify that the temperature inversion is 
responsible for the line doubling, the same line was computed with alt 
velocities set to zero. The result is shown in Figure V-ll. The central 
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reversal is characteristic of a line core formed in a region in which 
the source function increases outward. While a more detailed, non-LTE 
calculation would probably not show as large a central reversal since 
the source function in the line core could be smaller than the Planck 
function, the line doubling could still occur for the strongest lines. 

In fact, Thomas suggested that a temperature inversion was responsible 
for the "Cheshire Cat" lines (see Kraft 1967, pg. 239). The asymmetry 
of the line core can be caused by purely geometric effects as illus- 
trated in Figure V-12. This profile was computed with a constant 
velocity of expansion of 10 km s This profile is nearly identical 
to that in Figure V-10 except for the direction of the asymmetry caused 
by the different signs of the velocities used in the two cases. 

Xt is now possible to explain why the "Cheshire Cat" lines in 
and Call K are displaced more from the primary component than those 
of the metal lines. Since the temperature inversions in RDT occur 
only above t = 0-01, only the part of the line formed above this point 
in the atmosphere can show the central reversal. A large part of H 

CL 

is formed in this region, and, therefore, the velocity difference 
between the two components will be greater than that of the metal lines. 

Although the models are inadequate for computing profiles of 

strong lines, H profiles were computed for a qualitative comparison 

with observations. No attempt was made to accurately reproduce the 

observed line shapes. In fact, the H lines shown in Figure V-13 

CL 

were calculated using a Voigt profile instead of the more accurate 
Stark broadening. The lines have flat bottoms because the models do 
not extend to small enough optical depths and because of the LTE 
source function used. Qualitatively, though, the agreement with the 



profiles published by Grenfell and Wallers tein (1969) , Wallerstein 
(1972), Rodgers and Bell (1968), and Schmidt (1970) is excellent. 

While some of the irregularities in the profiles are associated with 
temperature fluctuations caused by zoning effects, several of the 
profiles indicate the presence of temperature inversions. In 
particular, the profile at $' ® 0.55 shows the effect of the shock 
wave that produces the second bump on the light and velocity curves. 

One of the problems associated with "Cheshire Cat" lines is the 
threshold effect (Skalafuris 1974). In Cepheids the line splitting 
occurs only for the strongest lines, while in W Vir and RR Lyr stars 
the line doubling often occurs in weaker lines. If the splitting is 
caused by temperature inversions, two phenomena contribute to the 
threshold effect. First, as a pressure wave moves outward, it 
encounters a decreasing density. Conservation of mass requires that 
the wave must accelerate. Thus, the strongest shocks occur highest 
in the atmosphere. The energy dissipation at the shock front produces 
a temperature rise. Second, as shown in Chapter II, the time it takes 
an element of gas to return to radiative equilibrium following a per- 
turbation increases as the optical depth decreases. In the region 
where the weak lines are formed, the gas requires only a few seconds 
to return to equilibrium, while higher in the atmosphere where the 
cores of the strong lines are formed, the time it takes the gas to 
return to equilibrium can be several thousand seconds (Whitney 1967). 

In- this time the shock will have travelled about half a pressure scale 
height and crossed several zones in the model producing a reasonably 
thick layer with an elevated temperature. The splitting of weak lines 



in population II variables is, then, indicative of shocks in their 
photospheres. 

ThQ structure of the shock front will be modified 

by the use of an artificial viscosity pressure, though. Instead of 
having a temperature discontinuity at the shock front, the temperature 
rise will be spread out over 3 or 4 zones. However, ignoring the 
effects of radiation from the heated gas at the shock front, the 
total temperature rise across the shock should be the same whether 
or not artificial viscosity is used (Richtmeyer and Morton 1967). 

If radiation effects are included, the error in the temperature rise 
depends on the ratio of two characteristic times. If the time it takes 
the gas to return to thermal equilibrium, t R , is shorter than the time 
it takes the shock to move a distance equal to its thickness, t , then 
using artificial viscosity will lead to an underestimate of the tem- 
perature rise. If, on the other hand, t > t , the temperature rise 
across the shock front will be independent of the artificial viscosity. 

In the region of interest in RDT t ~ 5000 sec, and the thickness of 

K 

9 

the shock is 3 or 4 zones, AX * 5 x 10 cm. If the shock moves at 
roughly the sound speed of 10 km s , then t^ ~ t g . The temperature 
inversions in RDT are, therefore, underestimates. This error is 
unimportant, though, since the assumptions made in computing the 
atmospheric structure are not completely valid in the region where the 
shocks occur. 

Care must be taken when measuring velocities of strong lines. 

The practice of treating each component of the line as a distinct 
layer of gas produces erroneous velocity curves. Neither component 
represents the motion of the atmosphere. Even if the splitting is 



not observed and the line core appears to be symmetric, there may be 
incipient splitting that has been masked by macroturbulence, micro- 
turbulence, or instrumental broadening. Velocities of strong lines 
should, therefore, be measured at some point in the wings of the 
line. This point will be discussed further in the next section in 
which the velocity curves are discussed. 

D. Velocity curves . 

1. Ratio of pulsational to radial velocity . 

In order to interpret radial velocities obtained from measure- 
ments of spectral lines, the center to limb variations of the 
specific intensity must be taken into account. In effect, the 
center to limb variations reduce to a determination of p, the ratio 
of the pulsation velocity of the atmosphere, v^, to the measured 
radial velocity of the line, v^. The value of p differs from 1 
because the observed flux profile is a weighted sum of intensity 
profiles with different line of sight velocities. (See Figure V-6). 

The early attempts at calculating p (Shapley and Nicholson 
1919; Carroll 1928; Getting 1935; van Hoof and Deurinck 1952) all 
made the assumption, either implicitly or explicitly, that the 
Doppler shift of the line was less than the line's Doppler width. 

This assumption is not always valid in Cepheids since the lines are 
typically 5 to 10 km s ^ wide while the pulsation amplitude is of 
the order of 20 to 30 km s In effect, the assumption reduces to 

the statement that the minimum of the flux profile is simply the 
weighted average of the minima of the intensity profiles from p. = 1 
to \i = 0. In the case shown in Figure V-7, it is clear that intensity 
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profiles from jj, < 0.4 have little effect on the location of the 
flux minimum. If the limb darkening law is given by 

ifei , (i-p) + pp, (v-1) 

KD 

this assumption leads to 

P . tzJt (V-2) 

4 - P 

Since equation (V-2) gives ff s P 5 ff for 0 * several 

investigations have concluded that the calibration of p has only 
& small effect on the computed velocities. Since the basic 
assumption used to calculate equation (V-2) is not valid for 
Cepheids, Parsons (1972) recalibrated p and showed that p is a 
function of y, the ratio of the velocity of the atmosphere, v pS to 
the observed width of the line including instrumental broadening v^. 

Since RDT differs from the models used by Parsons, a similar calcula- 
tion was performed. 

Three models were selected for the study, STB, $ ** 0.2, and 
$ s= 0.7, corresponding to mean, minimum, and maximum radius, 
respectively. For each model, p was computed for a set of veloci- 
ties and values of log gf. The effect of instrumental broadening 
was included by convolving each profile with a Gaussian slit 
function 0.041 wide, while rotation and macroturbulence were in- 
cluded by convolving the profiles with the rotation broadening 
function given by Huang and Struve (1960) with v sin i 8:1 10 km s . 
Although Kraft (1966) has concluded that Cepheids probably do not 
rotate, Abt (1958) has shown that this rotation cannot be distinguished 



from a macroturbulence of 7 km s t a value characteristic of Cepheids. 
The bisector of each profile was constructed and used to determine the 
velocity at the minimum, half intensity, and 1/e points of the profile, 
A total of 864 values of p was determined having an average of 
p - 1,31 ± 0.03. The value of p was found to be insensitive to the 
model, changing by only 0.02 between the models at minimum and maximum 
radius. If only those cases with line widths of 7 km s ^ are consid- 
ered, the average p = 1,30 ± 0.02, Since most Cepheids have line 
widths of this order, a constant value, p » 1.30, was adopted. The 
error in p introduces only a 2 % error in the velocity. 

Although a constant value of p was adopted. Figure V-14 shows 
some interesting correlations. The abscissae are y * v/v^ and the 
ordinates p. In all four graphs, the dashed lines are given by 
p = 1.37 - 0.03y taken from Parsons (1972), Figures (a) and (b) 
show p(y) for a weak and strong line, respectively, measured at the 
minimum of the profile while (c) and (d) represent the same lines 
measured at the half intensity points. In both cases, the strong 
line has a weaker dependence on y than the weak line. In addition, 
when the velocity is measured at the half intensity point, p is nearly 
independent of y, The near constancy of p when the velocities are 
measured at the half intensity point combined with the problem of 
incipient splitting of the cores of the strong lines indicates that 
the half intensity point should be used to determine velocities. In 
particular, the practice of using the line core when it is symmetric 
and the wings when the core is distorted will produce spurious scatter 
in the velocity curves. 



2» Comparison of hydrodynamic and line profile velocity curves . 

If the line core cannot be used to find the velocity of the upper 
atmosphere, the problem of detecting velocity gradients becomes more 
difficult. Figure V-15 shows how a velocity characteristic of the 
line core can be found. The dashed lines are extensions of the wings 
from points near the line core. The velocity determined from the 
point of intersection of the extrapolated wings agrees with the 
velocity of the upper layers of the model with an accuracy of about 
10X. In the example shown with p = 1.30, the predicted value is 
14.5 km s” 1 , while the model gives 13 km s * for the region in which 
the line core is formed. 

The easiest way to compare the velocities of the models with the 
velocities determined from the profiles is to compare velocity curves. 
Figure V-16 shows the velocity curves for a weak, intermediate, and 
strong line represented by a dashed, solid, and dot-dashed line, 
respectively. To make the comparison as representative as possible 
these velocities were measured at the minimum of the unconvolved 
profiles. The arrows indicate phases at which the strong line has a 
split core. Only the deeper component was measured leading to the 
deviations from the other velocity curves. The open circles show 
the velocity of H a determined by the extrapolation of the line wings 
discussed above. A comparison with Figure III-7 indicates that 
measurements of the line profiles underestimate the velocity differ- 
ences between the optically thin and optically thick layers. This 
result is not surprising since the line profiles represent averages 
over a range of continuum optical depths. Note, however, that the 

phase shift between K and the metal lines near $ = 0,5 is nearly the 

a 
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same as the phase Shift between t « 10 and T - 0.1 at the same 
phase. In addition, the feature on the velocity curve near $ = 0.2 
can be seen only for the strongest lines in agreement with Figure III -7, 
3, Center of mass velocity . 

One of the problems facing an observer studying Cepheids is the 
determination of the center of mass velocity of the star, v^. This 
value is needed both for kinematic studies (Kraft and Schmidt 1963; 
Wielen 1974) and for determining the radius of the star using the 
Wesselink (1946) method. Oke, Giver, and Searle (1962) have estimated 
that an error of 1 km s ^ in results in an error of 10% in the 
derived radius. 

The center of mass velocity is normally defined by finding 
Vq such that 


/Vv o )dt=0 . 

0 


<V-3) 


Due to changes in the opacity scale during the pulsation, the velocity 
curves do not follow a given element of gas. There is no guarantee, 
therefore, that a strict application of Equation (V-3) will give an 
accurate value of v^. 

In order to check the importance of this effect, velocity curves 
for a number of cases were constructed, and Equation (V-3) was used 
to find Vq. Table V-3 summarizes the results. In this table, MIN, 

1/e refer to the part of the profile used to find the velocity; 

N, G, R, B refer to the unconvolved profile, convolved with a 
Gaussian slit 0.04A wide, convolved with a rotation broadening 
function corresponding to v sin i = 10 km s \ and convolved with 
both broadening functions, respectively. The weakest lines were 
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completely washed out by the rotation broadening function and these 
cases are omitted from the table. The velocity curves determined 
from the minima of the strong lines have also been omitted due to 
splitting of the line core. It can be seen that v Q can be found to 
an accuracy of better than ± 0.4 km s by using Equation (V~3) . 

While this error will introduce a 4% error in the radius determination, 
it is rare that observations can be made to this accuracy. 

Inspection of these velocity curves reveals only small differ” 
ences among them. The adopted curve. Figure V-17, represents a 
typical velocity curve. The line used is moderately strong, varying 
from an equivalent width W * 88 n& to W = 153 mA. A macroturbulence of 
7 km s"^ has been assumed and the resultant profile was convolved with 
a Gaussian slit corresponding to 2 i / nan. In view of the discussion 
earlier in this chapter, the velocity was measured at the half inten- 
sity point on the profile. This velocity curve has been adopted as 
representative of the observed velocity curves of the metal lines 
and will be the only one used in the next chapter. 
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Table V-l 

Radial velocity in km s * observed 
from computed profiles 


log a 
a(km s 1 ) 

log m 0 

-1 

5 

-1 

10 

-2 

5 

-2 

10 

-3 

5 

-3 

10 

-1.000 

2.9 

5.4 

3.3 

5.4 

2.1 

5.4 

0.204 

5.4 

6.6 

5.4 

9.3 

4.2 

9.3 

1.408 

6.7 

13.3 

6.6 

14.7 

6.6 

14.7 

2.612 

10.8 

14.7 

9.3 

14.7 

6.6 

14.7 
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Table V-2 


Micro turbulent velocities in 
hydrodynamic model 


$ 

A log W 

*D 

v §/ v 

/V STB 

^ V STB 

1 

A V 

.00 

.12 

1.32 

0.97 

0.88 

1.39 

3.6 

.05 

.10 

1.26 

0,99 

0.77 

1.22 

1.6 

.10 

,06 

1.15 

1.01 

0.56 

0.89 

1.2 

.15 

.20 

1.58 

1.02 

1.22 

1.93 

5.8 

.20 

.16 

1.45 

1.05 

1.03 

1.63 

4.7 

.25 

.08 

1.20 

1.05 

0.62 

0.98 

0.8 

.30 

.14 

1.38 

1.05 

0.92 

1.46 

2.1 

.35 

.14 

1.38 

1.06 

0.92 

1.46 

2.0 

.40 

.10 

1.26 

1.05 

0.73 

1.16 

1.0 

.45 

.12 

1.32 

1.04 

0.84 

1.33 

2.4 

.50 

.18 

1.51 

i.02 

1.12 

1.77 

3.8 

.55 

.08 

1.20 

1.01 

0.66 

1,04 

1.6 

.60 

.20 

1.58 

1.01 

1.22 

1.93 

3.7 

.65 

.08 

1.20 

0.99 

0.67 

1.06 

1.6 

.70 

.04 

1.10 

0.98 

0.48 

0.76 

0.7 

.75 

.06 

1.15 

0.97 

0.59 

0.93 

1.0 

.80 

.08 

1.20 

0.97 

0.69 

1.09 

1.4 

.85 

.08 

1.20 

0.96 

0.69 

1.09 

2.3 

.90 

.06 

1.15 

0.96 

0,60 

0.95 

0.4 

.95 

00 

o 

• 

1.20 

0.96 

0,69 

1.09 

2.3 





Table V -3 


Error In center of mass velocity (km sec "*) 


MIN h ' 1 /e 


log gf 

N 

‘ c 

R 

B 

N 

C 

R 

B 

N 

C 

R 

B 

-4 

0:2 

0.4 

— 

— 

0.4 

0.4 

— 

— 

0.2 

0.2 

— 

— 

•3 

0.0 

0.3 

— 

— 

0.4 

0.3 

— 

— 

0.4 

0.3 

— 

— 

-2 

0.1 

0.1 

0.0 

0.0 

0.1 

0.1 

0.1 

- 0.2 

0.1 

0.0 

0.0 

0,0 

-1 

0.4 

0.0 

- 0.2 

- 0.2 

- 0.2 

- 0.2 

- 0.2 

- 0.2 

- 0.3 

■<• 0.2 • 

- 0.2 

- 0.2 

0 



— 


0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

0,1 

0.1 

1 





0,1 

0.1 

0.1 

0.1 

0.1 

0.0 

0.0 

0.0 


8 TT 
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Figure v-1. Curves of growth for log a * -1. — • • — — r 

Z = 0(C CQ ); , C - 5 km s" 1 (C^ 5 )? 

o o, a = 5 km s -1 <C a & ) ? • • , 

a * 10 km s” 1 (C aiQ ) ■ ^ is the microturbulent 

velocity and a is the velocity gradient parameter. 
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LOG a = -2 



LOG (ij/i) 0 ) 


Figure V-2. 


Curves of growth for log a - -2. Notation is the 
same as in Figure 1* 
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-LOG (W/X) 



LOG (tj/ ijo) 


Figure V-3. Curves of growth for log a = -3. Notation is the same 
as in Figure 1. 
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Figure V-5. 


Micrpturbulent velocity, vs. a measure of the velocity 
gradient, AV. The least squares line is shown. 
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Figure V-6. Schematic illustrating effect of angle integration on 
a line profile in an expanding atmosphere. 
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Figure V-7. Intensity profiles of lines in an expanding atmosphere. 

The curves are labelled with the corresponding value of 
p. •« cos 9 . 
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Figure V-8. a) Profile of a weak line in a moving atmosphere without 
velocity gradient. The dashed line is the bisector. 

b) Same as (a) for a strong line. 
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X(A) 


8 V-9. Profile of a strong line in an atmosphere with a velocity 

g dient. The dashed line is a straight line drawn 
vertically from the minimum of the profile. The bisector 
Of the profile is also shown. 
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Figure V-10. Line profile showing a "Cheshire Cat" line. 
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Figure V-12. Same case as Figure V-10 but with a constant velocity 
in the atmosphere. 
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Figure V-I4. Ratio of pulsation to observed radial velocity, p, vs. 

ratio of pulsation velocity to half width of line, y . 
Dashed line is Parsons 1 (1972) relation. 

a) Weak line measured at minimum of profile 

b) Strong line measured at minimum of profile 

c) Weak line measured at half intensity point 

d) Strong line measured at half intensity point 
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Figure V- 15'.' Method for computing velocity characteristic of line core 

when core is spLit. Velocity is measured from intersection 
of dashed lines. 
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Figure V-16. Velocity curves. intermediate strength line, weak 

line, . strong line, o H determined using method 


described in text. 



Figure V-17. Adopted velocity curve measured from line profiles 
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CHAPTER VI 

DETERMINATION OF CEPHEID RADII 

A. Radius determination methods . 

X. Bolometric radius . 

There are two ways to find the radius of a Cepheid from 
observed quantities. The first method involves solving 

t = 4nR 2 <JT* ff . (VI-1) 

for the radius. The radius found in this way will be referred to 
as the bolometric radius. The mean luminosity of the star is 
found from the period-luminosity law while, in the simplest appli- 
cation, the temperature is calculated from the color. There are 
obvious drawbacks. A change of the zero point of the period- 
luminosity law of rffl changes log R by 0.02, about 47o. In addition, 
the solution is very sensitive to errors in which, if broad 

band colors are used to define T depends on the assumed reddening. 
Whitney (1955) attempted to improve this approach by using model 
atmospheres to determine the flux of the star and finding log R*/Rq 
from 

Ife - M* * 2.5 log 

where is the sensitivity function of the filter used to find 
Oke (1961 a, b) made a further improvement by making absolute flux 
measurements in 50 A bands. He then compared these fluxes to those 
computed from model atmospheres to find T While both of these 

approaches remove some of the errors inherent in the color-T ^ 


7 5 — l + 5 log (VV (VI -2) 

I fx\» 

l j 
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calibration, they require knowledge of the line blocking. An 
error of only 50°K in T ff results in an 0.01 change in log R. In 
addition to these systematic errors, the gravity variations, which 
can be as large as a factor of 10 (Parsons 1971 a, b; see also 
Chapter IV) will produce an error that varies during the pulsation 
cycle. 

2. Baade and Wesselink radii . 

One way to find the radius that avoids these problems is the 
Wesselink (1946) method. Based on an idea of Baade (1926), the 
method uses both the light and velocity curves to find the radius. 
Baade proposed using the change in brightness to give the ratio of 
the radii from 

- M 2 * 5 log R 2 /R x + 10 log T 2 /T^ , (VI-3) 

and the velocity curve to give 



Baade's method requires knowledge of the color temperature law. It 
was not surprising when Bottlinger's (1928) attempt to find the 
radius of £ Gem failed since he assumed the star radiated like a 
black body. Becker (1940) improved the situation by assuming only 
that a single valued color-T^^ law existed and obtained radii for a 
number of Cepheids. 

Wesselink (1946) removed the problem of the color-T^^ law when 
he suggested choosing two phases at which the star has the same 
color. If it is assumed that equal color implies equal temperature, 
the last term on the right of equation- (Vl-3) vanishes, and there is 



no need to know even the form of the color-*T -- relation. The 

eft 

Wesselink method also has the advantage that it is independent of 
both the zero point of the period- luminosity law and the interstellar 
reddening . 

There are two other assumptions inherent in the Wesselink method. 
First, the mass depth of the line forming region is assumed constant 
with phase so that the observed velocity curve follows a given element 
of gas. Second, it is assumed that the ratio of the radii of the line 
forming and continuum forming regions is constant with phase. These 
assumptions will be examined below. 

There are, of course, problems with the method. Accurate velocity 
curves are required and any change in p, the ratio of pulsational to 
observed radial velocity, produces a systematic error in the radius. 

In addition, loops in the (U-B) - (B-V) diagram indicate that equal 
color does not necessarily imply equal temperature. It is also known 
that the opacity scale changes during the cycle so that different 
elements of gas are observed at different phases. Velocity gradients 
measured in Cepheid atmospheres by Sanford (1956) and Dawe (1969), 
among others, indicate that the ratio of radii of the photosphere 
and reversing layer changes during the pulsation cycle but the size 
of this change is not known. Fernie and Hube (1967) have also shown 
that small errors in reducing the light and velocity curves to the 
same epoch produce large errors in the derived radius. In spite of 
these difficulties excellent results have been obtained. 

B. Calculated radius of the hydrodynamic model . 

1, Wesselink radius . 

Most papers reporting Cepheid observations include a section 
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computing the radius by Wesselink's method. In accordance with this 
tradition, the results of the previous chapters can be treated as 
observations and the radius of RDT calculated using Wesselink's method. 
The light and color curves of Chapter XV and the velocity curve of 
Chapter V have been adopted as if they were observations of a star. 

As pointed out previously, zoning effects limit the accuracy of the 
light and color curves to about oTftf which translates into an 6% 
random error in the radius. The integral of the adopted velocity 
curve indicates an error of 0.2 km s” 1 in the center of mass velocity 
of the star. Combined with a random error of 2% from changes in p, 
the minimum error expected is 10%, 8% random and 27* systematic. Since 
Fernie (1968) was dealing with less accurate velocity data than that 
used here, hie adopted error of 107* appears to be an underestimate. 

Figure VI- 1 shows the radius curve derived from the adopted 
velocity curve (dashed line) and the radius of T * 1 taken directly 
from the models (solid line). The amplitude of the radius curve 
determined from the lines is about 77* larger than the variation of 
the photospheric radius. This error will result in an overestimate 
of the radius amplitude but should have only a small effect on the 
computed value of the mean radius. It appears, therefore, that errors 
introduced into the mean radius by changes of the opacity scale with 
phase and by velocity gradients in the atmosphere are small. 

Figure VI-2 shows the results of the Wesselink calculation 
performed by taking pairs of points with equal (B-V), (V-R) , and 
(R-I). The solid line was taken directly from the models. The 
error bar shows that the expected error of a single measurement is 
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comparable to the total radius variation. Therefore, all that can 
be derived from these values is the mean radius. STB has a radius 
of 71.7 Rq while the Wesselink calculations give 58.5 ± 4.0, 73.5 ± 
8.7, and 67.3 db 5.2 for (B-V) , (V-R), and (R-I) , respectively, 
where the quoted errors are the standard deviations computed from 
several radius determinations. The agreement is satisfactory except 
for the (B-V) curve which is most sensitive to changes in It 

thus appears that the Wesselink method can be used to find the mean 
radius of a Cepheid with an accuracy of about 107« only if very high 
accuracy observations are available. 

2. Baade radius . 

At the time Wesselink published his modification of Baade's 
method, the color-T^^ relation was not known accurately. Since 
then the relation has been calibrated for several different color 
systems. In Chapter IV, it was shown that the color- T ef f relations 
derived for RDT agree with those derived from observations. Figure 
VI -3 shows the values of R/Rq, where is the radius at $ = 0, 
derived using Baade's method. The notation is the same as in Figure 
VI-2. Using Baade's method introduces an additional error. A change 
of 0?02 in the color results in an error of 2% in the radius. Baade's 
method does not depend on the zero point of the color-T^^ law or 
that of the reddening law. Only the slopes of these relations enter 
the calculations. 

The agreement with the model is quite good except near maximum 
light ($ - 0.35). The radii derived from (R-I) are more accurate 
than those of (B-V) and (V-R) , especially near minimum light. Com- 
bined with Figure VI-1, mean radii in solar units of 77.2 ± 7.8 from 


(B-V), 79.6 ± 4.7 from (V-R), and 73.6 ± 5.5 from (R-I) are derived 
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If the points near maximum light are removed, these values become 
71,0 ± 5.5 (B-V) s 74.2 i 3.8 (V-R), and 73.6 ± 2.8 (R-I)o Again the 
agreement is satisfactory, and there appears to be no need to select 
only phases of equal color unless the reddening is not known. 

3. Bolometric radius . 

One further set of radius determination was made, this time 
using Equation VI-1. Temperatures were computed from the color-T e £^ 
relations derived in Chapter IV * The bolometric magnitude of the 
sun, = 4.72, and the solar effective temperature, = 5800°K 9 

were taken from Allen (1963). In this case an error of 0.01 in the 
zero point of the color-T^^ relation introduces a systematic error 
of 57« into the radius determination. A random error of 8% and a 
systematic error of TL can, therefore, be expected. The mean radius 
derived is 69 ± 2 regardless of which color was used. The agree- 
ment is again satisfactory, and this method can be used for stars 
that do not have accurately measured velocity curves. The radius 
determinations are summarized in Table VI-1. 


C. Method of Woo ley and Savage . 

Some comments are in order on the modification of the Wesselink 
method proposed by Wooley and Savage (1971) and Wooley and Carter 
(1973). Their main point is that the color depends on g & ££ as well 
as T e ^ f . Thus, it is necessary to include the luminosity of the star 
in the calculation, which means both the mass and radius can then be 
found. They also use a velocity function instead of the observed 
velocities which are often too poorly determined to be useful. 

Figure VI -4 compares the adopted velocity curve of RDT (solid line) 



with the velocity function of Wooley and Carter (1973) scaled to a 
semi-amplitude of 15 km s * (dashed line). There is a reasonably 

good agreement between the two curves except near the second bump. 

Wooley and Savage make two unnecessary assumptions, however. 

First, they make the assumption that A R/R is small enough that 

In (1 + Ar/r) « AR/R. Fernie and Hube (1967) have shown that this 

assumption leads to an underestimate of the radius of 5 to 10%. In 

particular, the radius of RDT would be underestimated by 11%. Second, 

they assume that the accelerations of the atmosphere are always small 

2 

compared to g s G M/R when in fact they are often much larger than 
g (Parsons 1971; see also Chapter IV). These two underestimates are 
partially cancelled by using p = 24/17 = 1.41, instead of the smaller 
value suggested by Parsons (1972) and in Chapter V, 

As shown in this chapter, though, the gravity variations have 
only a small effect on the derived mean radius, especially if (R-I) 
is used. Even the change of the opacity scale with phase produces 
only small errors in the mean radius if sufficiently weak lines are 
used to define the velocity curve. With the best observations it 
should be possible to achieve an accuracy of 10%, with any of the 
methods discussed above. In practice, however, the errors may be 
closer to 15 or 20% due to additional errors introduced by the lower 
accuracy of the observed velocities. 



Table VI-1 


Radius determinations of hydrodynamic model 
in solar unite 


Method 

B-V 

R ± s ,d. 

r stb " 71,7 

V-R 

R ± s.d. 

R-I 
R ± s. 

d. 

expected 

random 

error 

systematic 

Wesselink (1946) 

58.5 ± 4,0 

73.5 ± 8.7 

67.3 ± 

5.2 

8% 

2% 

Baade (1926) 

77.2 ± 7.8 

79.6 ± 4.7 

73.6 ± 

5.5 

10% 

2% 

all points 
some points* 

71.0 ± 5.5 

74.2 ± 3.8 

73.6 ± 

2.8 

10% 

2% 

2 4 

L « 4 TTR o T 

69.4 i 2.0 

69.5 ± 2,2 

69,2 dt 

1.6 

8% 

7% 


^excluding points near maximum light 
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Figure VI- 1. Change of radius, R-R Q , vs. phase from models, from 

adopted velocity curve. 
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line from models? 0 from (B-V ) , + from 
(V-R) , x from (R-I) . 
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Figure VI-4. Adopted velocity curve ( ) compared to velocity 

function of Woolley and Carter (1973) normalized to a 
semi-amplitude of 15 km s”l ( ) . 




CHAPTER VII 


SUMMARY, CONCLUSIONS, AND FUTURE WORK 
A. Summary and Conclusions . 

Hydrodynamic models of the atmosphere of a 12 Cepheid have 
been computed including the effects of radiative transfer in the 
optically thin zones* A new method of including radiative transfer 
in a standard Henyey type hydrodynamic code has been developed. The 
differences between using the diffusion approximation and the solu- 
tion of the transfer equation have been shown to be negligible ex- 
-2 

cept above T = 10 where temperature inversions occur only in the 
radiative transfer case. A study of the envelope of the full 
amplitude model indicates that the phase lag between maximum light 
and minimum radius increases continuously between the Hell and H 
ionization zones. The asymmetry of the light curve appears to origi- 
nate in the H ionization region. 

The Hertzsprung sequence has been examined and a mechanism 
presented to explain the occurrence of two bumps on Cepheid light 
curves. The bump occurring on the falling branch of the light 
curves of 7^ to 10^ Cepheids and on the rising branch of 10^ to 15^ 
Cepheids is caused by the Christy mechanism, i.e., a pressure wave 
which propagates into the star, reflects from the stellar core, and 
appears on the next pulsation cycle. The other bump which appears 
on the falling branch of 10^ to 15^ Cepheids is due to an atmospheric 
oscillation. This bump occurs when the natural pulsation mode of 
the atmosphere, which has a shorter period than the envelope, has a 
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large enough amplitude to generate a shock wave by comp res sing the 
hydrogen ionization region. When this shock reaches the surface 
it appears as a bump on the falling branch of the light curve. Since 
two mechanisms produce bumps on Cepheid light curves , masses derived 
from the phase of the bump may be unreliable* While it has been 
shown that this mechanism is consistent with the Herzsprung sequence, 
more models must be computed to fully study this atmospheric oscilla- 
tion mechanism* 

The hydrodynamic models were then used to compute UBVRI colors 
by treating the model at each time step as a snapshot of the atmo- 
spheric structure. A P-V diagram of the atmosphere constructed from 
the effective temperature and gravity at each phase indicates the 
destabilizing influence of the hydrogen ionization region. These 
effective temperatures and gravities were then used to compute line 
blocking coefficients which were applied to the monochromatic fluxes 
to give the colors of the models, including the effects of spectral 
lines. The color curves were smoothed to minimize the zoning effects, 
and it was shown that they reproduce the observed variation of light 
amplitude and phase of light maximum with effective wavelength. The 
color-T e ^ relations were computed and were shown to agree with 
those derived independently. It was then shown that the loops in 
the (U-B)-(B-V) diagram are most likely caused by a non thermal depen- 
dence of the continuous opacity. Mean colors of the model were computed 

if 

using three averaging schemes, and it was found that the intensity means 
of the magnitudes, (b)^ - (v)^, best represent the colors of the 
equilibrium model. The location of the equilibrium model in the 
H-R diagram indicates that the zero point of the Sandage and Tammarm 
period- luminosity, relation is too high by 0?2. 



Line profiles were then computed using the moving atmospheres 
from the hydrodynamic models „ It was shown that, although the 
velocity gradients in the atmosphere are not responsible for the 
observed microturbulence or variation of microturbulence, they 
can be used to explain the occurrence of supersonic micro turbulence. 
The total observed microturbulence was shown to be consistent with 
the linear sum of the classical microturbulence and that caused by 
the velocity gradients. 

The splitting of the cores of strong lines was shown to be due 
to shock induced temperature inversions in the Cepheid atmosphere. 

This mechanism explains why the splitting is observed only in strong 
lines in classical Cepheids and why the splitting in H and Call H 
and K is greater than in the strong metal lines. The splitting of 
the line core makes velocity measurements from the minimum of the 
profile unreliable, but a method for determining the velocity of 
the upper atmosphere was presented. The center to limb variations of 
the profiles were then studied, and it was shown that the commonly 
used ratio of the pulsation to observed radial velocity, p = 24/17, 
is too high. It was found that velocities obtained from the lines 
underestimate the velocity gradients present in the model atmospheres. 
It was also shown that the integral of the velocity over phase can be 
used to find the velocity of the center of mass of the star to the 
accuracy of observations in spite of changes in the continuous opacity 
scale during the pulsation. 

The adopted light, color, and velocity curves were then used 
to study various methods for determining the mean radius of a 
Cepheid. The Wes se link method was found to give radii accurate to 



©bout 10% if errors in the observed velocity curves are minimized . 

Baade r s method also produces accurate radius determinations if the 

currently accepted color-T eff relations are used. It was also shown 

2 4 

that the bolometric radius found from L * 4tr R o r is reasonably 
accurate and can be used for stars for which the velocity curve is not 
accurately known. However s the bolometric radius is susceptible to 
systematic errors introduced by errors in the reddening, the zero point 
and slope of the color-T gff law> and the zero point of the period- 
luminosity relation* The Wesselink method is independent of these 
systematic errors, and Baade* s method is affected only by errors in 
the slope of the color- T ^ relation. All these errors can be reduced 
by vising R-I instead of B-V. 

B. Future work . 

There are several modifications to the models that would improve 
the agreement with observations, i.e. s nongray radiative transfer, proper 
treatment of line blanketing, convective energy transport. The most 
important of these changes is the inclusion, of convection. While it is 
true that convection can carry only a small part of the flux due to the 
low densities of Cepheid envelopes, the destablizing influence of the 
hydrogen ionization region is sensitive to changes in the temperature 
gradient. The atmospheric pulsation modes are also dependent on the 
temperature structure of the hydrogen ionization region. Since it takes 
a ’convective element about 0.1 period to move one pressure scale height, 
theory of time dependent convection is needed. 

Some of the assumptions made in calculating the models need to be 
investigated. One of the most common of these assumptions is that the 



lower boundary of the model can be kept fixed. Both the velocity 
and kinetic energy asymptotically approach zero in deep envelope 
models, but, as shown in Figure VII- 1, a 3-D plot of momentum vs. 
mass point and phase, the momentum does not. The atmosphere has 
very little momentum due to the low density, but the momentum is 
being arbitrarily forced to zero at the base of the envelope by 
the zero velocity boundary condition. While it is not clear that 
this constraint is significant, a model should be computed with a 
free lower boundary to see if this assumption affects the observables 
Another interesting problem is the inhomogeneity of Cepheids. 
Examination of a catalog of light curves indicates that light curves 
of Cepheids with nearly the same period show striking differences. 

For example, £ Gem with a period of 10^15 has a low amplitude, nearly 
sinusoidal light curve, while £ Dor with a period of 9?84 has a 
large amplitude, asymmetric light curve with two distinct bumps. Our 
understanding of these differences is limited by the common practice 
of computing a model and finding a Cepheid with a similar light curve 
More would be learned about these differences if a grid of models 
were computed to fit a single star as is done for stellar atmospheres 
Such a grid could also be used to study the origin of the bumps on 
the light curves. 

The hydrodynamic model atmospheres can also be improved by 

including more optically thin zones extending to smaller optical 

depth. Ideally these models should extend to a Rosseland mean 

-6 

optical depth of 10 to allow the strongest lines and the possible 
formation of a chromosphere to be studied. This extension of the 



model will probably require that the plane parallel assumption 
be dropped. 

The line profile calculations can also be improved by including 
non-LTE effects, especially if "Cheshire Cat" lines are to be studied. 
The ratio of pulsation to radial velocity, p, should be studied 
further since its variation with line width and line strength is not 
understood. In particular, the variation of p with pulsation velocity 
should be included when constructing velocity curves as demonstrated 
ly Duquesne and Schatzman (1955) » Until the variation of p is under- 
stood errors in the radii determined by the Wesselink or Baade methods 
cannot be reduced much below 10% and the period -radius relation cannot 
be used to detect stars pulsating in overtone modes. 






















APPENDIX A 


COEFFICIENTS OF THE INCREMENTS— DIFFUSION APPROXIMATION 

The nonvanishing coefficients of the increments as derived by 
G. S. Kutter and W. M. Sparks used in the diffusion approximation cal- 
culations are listed. These coefficients are denoted by the symbol DET 
followed by two numbers. The first of these refers to the differential 
equation (II- 1 to II-5). The second number runs from 0 to 10, where 0 
refers to the inhomogeneous term, and 1 to 10 identify the coefficients 

of 8v i-l* 5B i-l* SR i-l’ 6W i-l/2’ 52 i-l/2 , 6v i’ 6B i* 5R i’ 5W i+l/2* 5Z i+l/2’ 
respectively. For instance, the conservation of mass equation (II-l) 

has the form 

DET 13 • 6R ±-1 + DET 14 * + DET18 - &K ± * DET10. 

i * 2: 

The coefficients are identical to those listed below for 1*3, ...,N 
except for 

DET13 * DET32 * 0. 

1*3, 

DET13 - COKl^/2 r[_ v 

DET14 - V 1-1/2 

DET18 * - co ^ 1 i „i/2 r i’ 

DET10 = ~V ± ~i/2 + I C0N1 i~l/2 (r i * 

C0Nl^ 1/2 * -4ir[M {) (Q i - Q^)] J 
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APPENDIX B 


COEFFICIENTS OF THE INCREMENTS— 

RADIATIVE TRANSFER 

The coefficients of the increments that are changed to include radiative 
transfer effects are listed. The notation is defined in Appendix A. 
i * 2, N: 


DET40 « B - B , 

1 \ 

DET44 * F1 ± (DET44) d , 

DET45 = Fl i (DET45) D 
DET48 » 0 , 

DET49 * F2 i (DET49) ]) , 

DET410 = F2 i (DET410) D , 

\ " 16,72 r eff / 1 1 viv ’ 


r 


eff 


r^ if £ 10 , 

r if t. < 10 , 

a i 


r 

a 



radius of deepest zone with t < 
1 if 10 , 

l-e~ ATl if t ± < 10 , 

7 1 if 10 , 

-At , 

1-e 1+1 if x < 10 , 

i 


10 


At. 


2^ T i ~ X i+2^ * 


* 
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i - H + Is 


DET40 * , 


DET44 - F2 n+1 (DET44) d , 


DET45 = F2 n+1 (DET45) d 


DET48 -* 0 , 


F2 N+1 


* 1 - e 


A subscript D denotes the quantity was computed using the formula in 


Appendix A. 



APPENDIX C 


ELECTRON PRESSURE ITERATION PROCEDURE 


The electron pressure P^ is important for determining the emergent 
spectrum since it affects the ionization balance through the Saha equation. 
P e> though, also depends on the number of free electrons, P^ = N^ kT. 
Therefore^ an iterative procedure must be used to find P g from the known 
temperature and density. The iteration used here can be written as 


a l + a 2 

[a 2 /P e (k " 1>] - 1 


where k counts the number of iterations. 


p (°) „ 1 p „ 1 QkT 

e 2 N 2 ym^ 


a 1 = P M V y.YjX. , 

1 N i l l 

i«l 

No 2 

» P„ £ Y.(Y. - X.Z.)/X. 

2 N ^ - l x x x' x 


y, = number fraction of element i 
x 


x. = i + E f . . 

1 j-1 13 


Y. = E jf . . 

i j = l J 13 


z. = ^ j 2 f _ 
x j = 1 J xj 




N i,n+1 


13 1 N. 

n=l x,n 


Convergence of this iteration is quadratic and global, 
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Plate 


Lu mi no si l ' 


mass point and 
areas , 


I a r g. e 1 u ra i n o s i t y ; d a r k 
envelope is at the top 
are shown. 


phase. Bright areas represent 
low luminosity* Base of the 


surface, at the bottom. Two periods 


^4 


Plate ! I . Velocity vs, mass point, and phase. Bright areas indicate 
expansion; dark areas contraction. Base of the envelope 
is at the top; surface at the bottom. Two periods are shown. 
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